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ABSTRACT 


Several  types  of  seismic  signal  detectors  have  been  tested  on  NORSAR 

and  Pinedale,  Wyoming  signals  buried  every  ten  minutes  in  approximately 

22  hours  of  phase-scrambled  noise.  All  plausible  detectors  including  human 

analysts  operate  within  0.1  m^  of  the  mean  value  including  the  present 

on-line  "IBM"  STA/LTA  detector,  a  "Z"  detector,  detectors  with  fixed  or 

optimum  filters,  a  "deflection"  detector,  the  MARS  detector,  a  detector 

which  cascades  a  recursive  digital  filter  with  Walsh  transforms,  and  human 

analysts  viewing  unfiltered  Calcomp  plots.  The  best  detector  was  a  log-Z 

2 

transform  aided  by  a  recursively  updated  optimum  S/N  filter.  Among  the  worst 
were  human  analysts  when  operating  in  the  "pick  everything"  mode,  the 
deflection  detector,  and  a  prediction-error  type  detector.  The  human  analyst 
could,  however,  improve  to  be  nearly  the  best  if  sufficient  care  was  taken. 
With  a  modest  amount  of  care,  any  detector  can  be  modified  to  perform  well 
as  an  element  of  an  operational  system  and  the  choice  of  single  channel 
detector  can  be  made  on  other  grounds  than  pure  detection  such  as  avail¬ 
ability  or  computation  speed. 

In  this  report  also  we  outline  the  features  of  the  existing  on-line 
detector  at  the  SDAC  which,  in  addition  to  pure  detection,  also  classifies 
events  as  spikes,  drop-outs,  regionals,  and  teleseisms  and  calculates  a 
refined  start  time. 
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INTRODUCTION 


As  discussed  in  more  detail  in  Appendix  I,  Freiberger  (1963)  derived  the 
optimum  single-channel  prefilter  and  detection  procedure  using  Neyman-Pearson 
arguments,  and  Vanderkulk  et  al  (1965)  made  use  of  this  theory  in  the  design 
of  the  operational  detector  developed  for  LASA.  Vanderkulk  et  al  evaluated 
by  means  of  theoretical  analyses  the  degradation  which  would  result  from 
various  approximations  to  the  optimum  procedure  and  concluded  that  the  dif¬ 
ference  would  be  small.  Further  research  on  detection  proceeded  to  questions 
of  array  detectors  which  might  be  resistant  to  false  alarms  due  to  spikes  and 
locals,  Blandford  (1970,  1972,  1974),  to  general  questions  concerning  array 
and  network  detection,  Wirth  et  al  (1971),  Shumway  (1971),  Ringdal  et  al  (1972), 
Blandford  and  Wirth  (1973)  and  Wirth  et  al  (1976)  ;  to  methods  of  holding  the 
false  alarm  rate  constant,  Lacoss  (1972),  and  to  various  ad-hoc  methods  of 
recognizing  spikes  and  local  events,  von  Seggern  (1977). 

Despite  Freiberger' s  result  a  number  of  workers  continued  to  program  and 
test  various  "intuitive"  detectors,  e.g.  Allen  (1978),  Shensa  (1977),  Savino 
(1979).  This,  perhaps,  derived  from  the  fact  that  there  has  never  been  a 
direct  comparison  on  the  same  data  base  of  the  optimum  detector  to  other 
single  channel  detectors.  To  remedy  this  situation  we  have,  as  discussed  in 
Appendices  II  and  III,  constructed  test  data  tapes  on  which  various  research 
workers  may  exercise  their  detectors  and  via  which  different  detectors  may  be 
compared. 

In  the  rest  of  this  study  we  first  outline  the  capabilities  of  the  com¬ 
puter  program  which  we  have  written  to  test  various  broad-band  detectors, 
and  then  compare  the  results  (on  application  to  the  test  tape)  of  STA/LTA 
detectors  with  various  averaging  times,  rectification  or  squaring,  band-pass, 
Weiner,  auto-regressive  or  prediction-error,  and  optimum  filters,  and  Z-trans- 
formation  with  or  without  the  log-normal  assumption  (Shensa,  1977).  The 
conclusion  from  all  this  is  that  the  standard  band-pass  filter  with  rectifica¬ 
tion  is  within  0.1  magnitude  unit  in  performance  of  the  optimum  process.  We 
also  discuss  the  performance  of  the  analyst  on  the  data  when  presented  with 
the  data  in  standard  unfiltered  develocorder  format  and  compare  with  results 
of  other  contractors'  work  on  these  tapes. 
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OUTLINE  OF  COMPUTER  PROGRAM  CAPABILITIES 

The  basic  mode  of  operation  of  the  360/44  research  program  NDP  is  that 
of  a  conventional  STA/LTA  detector.  The  data  may  be  read  in  either  from  a 
subset  tape  (INVSC=0)  or  from  the  VSC  detector  tapes  discussed  in  Appendices 
II  and  III,  (INVSC=I).  If  the  data  are  from  subset  then  other  parameters 
specifying  the  seismogram  number,  channel  and  the  noise  and  signal  windows 
are  needed,  but  we  shall  not  consider  this  case  further  in  this  study.  When 
the  data  are  read  in  they  may  then  either  be  decimated  (120=1)  or  not  (120=0); 
in  this  study  the  Pinedale  data  is  decimated  and  the  NORSAR  data,  originally 
at  10  sps,  is  not.  The  sampling  rate  of  the  resulting  data  is  specified,  in 
this  study  it  is  always  10  (SR=10) .  After  filtering,  which  we  will  discuss 
later,  the  filtered  data  stream  is  passed  to  the  main  detection  subroutine 
(STh/LTA).  The  filtered  data  are  either  squared  (IASQ=1)  or  the  absolute  value 
is  taken (IASQ=0) and  the  result  is  averaged  into  "YJ"  values  over  non-overlap¬ 
ping  windows  of  SA  points.  A  sliding  average  of  RP  of  these  values  is  computed 
and  constitutes  the  STA.  For  a  10  second  STA  the  parameters  used  are  SA=33  and 
RP=3  (thus,  actually  a  9.9  second  STA).  For  all  runs  in  this  study,  RP=3, 
and  SA  is  adjusted  to  give  the  appropriate  time  window.  The  LTA  is  recur¬ 
sively  updated  every  KR  times  an  STA  is  computed,  and  in  this  study,  KR=3, 
so  that  the  LTA  is  updated  every  time  a  fresh  STA  is  available;  therefore, 
at  time  intervals  which  increase  as  the  STA  window  length  increases.  This 
ensures  that  the  LTA  will  be  more  stable  than  the  STA.  The  LTA  is  updated 
with  a  recursive  filter  according  to  the  formula 

LTA  =  2SIG  *  STA  +  (1  -  2SIG)  LTA 

and  in  all  runs  in  this  paper,  SIG=-5.  After  detection  the  LTA  is  updated 
more  rapidly  with  SIGT  replacing  SIG  and  SIGT=-4  so  that  the  event  will 
"turn  off"  more  rapidly  to  ensure  that  later  phases  may  be  detected.  The 
fact  that  SIG=-5  and  RP=KR=3  implies  that  the  LTA  averages  over  a  window 
about  16-32  times  longer  than  the  STA  for  all  runs  in  this  study. 

The  threshold  for  detection  is  V  which  may  be  calculated  internally 
in  the  program  or  set  a  phsiofii  .  If  it  is  set  a  pviofu.  then  it  is  input  as 
VI;  and  the  turn-off  threshold  is  UI.  If  it  is  determined  and  updated 
internally  then  it  is  calculated  from  the  desired  false  alarm  rate  which  is 
specified  as  the  signal  window,  T  divided  by  the  desired  time  between  false 
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alarms  TI,  both  expressed  in  seconds.  The  actual  modes  of  threshold  calcula¬ 
tion  are  discussed  below. 

If  STA/LTA  >  V  for  QT  out  of  QP  successive  tests  then  a  detection  is 
declared.  The  on-line  SDAC  system  uses,  for  reasons  of  historical  continuity, 

QT=3  and  QP=3,  and  when  we  are  trying  to  exactly  duplicate  that  detector  or 
small  modifications  of  it,  these  parameters  are  used.  In  general,  however, 
in  this  paper  we  use  QT=1,  QP=3,  which  simply  amounts  to  detection  if  STA/LTA 
exceeds  V  even  once. 

Various  options  exist  for  filtering  of  the  raw  data.  There  is  provision 

for  recursive  filtering  and  several  sets  of  coefficients  for  bandpass  filters 

exist  in  the  program  in  data  statements.  Those  used  in  *:his  study  are  all  for 

10  sample  per  second  data  and  are  6-pole  (IRP=6)  Butterwurth  filters  for 

1.5-4  Hz  (suitable  for  NORSAR)  and  1-1.5  and  1-2  Hz  (suitable  for  Pinedale). 

Recursive  filtering  is  carried  out  by  subroutine  RFIL  if  IOPTF=0.  The  optimum 

filter  is  applied  if  I0PTF=1.  If  ISNEQ=0  the  signal  is  assumed  to  be  small 

2 

and  the  optimum  filter  is  S/N  ,  (see  Appendix  I).  If  ISNEQ=1,  then  S  is  ad¬ 
justed  in  amplitude  so  that  S/N  <_  1  but  is  just  equal  to  1.0  at  some  frequency. 
This  signal  level  should  be  just  detectable  and  is  therefore  the  signal  level 

for  which  we  optimize  the  detector.  Then,  this  S  is  used  to  compute  the 

2  2  1/2 

optimum  filter  S/[N(S  +  N  )  ].  Unless  otherwise  noted,  this  is  the  filter 

2  2  2 

used  in  this  study.  To  compute  Weiner  filters  S  /(S  +  N  ),  a  special  but 

simple  change  to  the  code  is  made  and  again  S  is  adjusted  so  that  S/N  £  1. 

The  signal  amplitude  spectrum  S  is  computed  as  the  product  of  a  flat  source 
spectrum  times  the  instrument  response  (The  instrument  response  is  designated 
in  a  data  statement  and  is  recovered  from  the  RSPCOM  disc  file.)  times  the 
function  exp(— ^  )  where  t*  is  set  to  be  0.3  unless  otherwise  noted  in  this  paper. 
A  flat  source  spectrum  for  P  waves  from  small  events  would  be  deduced  from 
almost  all  source  theories  and  observations  in  the  literature  although  for 
purposes  of  this  experiment,  in  which  signals  from  larger  events  were  buried 
in  noise,  a  corner  frequency  in  the  range  1-4  Hz  might  have  been  preferable. 

In  any  event,  we  have  adopted  the  flat  spectrum  as  the  simplest  possibility, 
and  the  one  which  would  be  correct  in  application.  Obviously,  higher  values 
of  t*  would  be  desirable  for  events  from  rift  zones  and  other  low-Q  source 
regions,  and  lower  values  would  be  suitable  for  events  from  shield  regions  to 
NORSAR.  We  shall  see,  however,  that  the  performance  is  not  overly  sensitive 
to  this  parameter. 
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The  noise  spectrum,  N,  can,  of  course,  be  determined  directly  from  the 
trace  under  consideration.  There  is,  naturally,  no  need  to  know  the  system 
response  to  calculate  the  noise  spectrum  since  the  noise  spectrum  as  we  see 
it  has  passed  through  the  system  response  already.  The  noise  spectrum  is 
obtained  from  the  first  IP  points  (IP=512)  of  the  10-minute  noise  plus 
signal  window,  and  is  smoothed  with  a  running  averaging  window  ISMU  points 
long  (ISMU=9).  If  ISL0G=1  then  the  log  spectrum  is  smoothed  but  in  this 
study,  ISL0G=0  and  the  amplitude  spectrum  is  smoothed.  As  we  shall  see  the 
opcimum  filter  determined  as  discussed  above  was  usually  slightly  outperform¬ 
ed  by  a  fixed  bandpass  filter  whose  cutoff  frequencies  were  determined  by 
examination  of  several  of  the  optimum  filter  shapes.  This  suggested  that  the 
optimum  filter  was  somewhat  unstable,  not  being  averaged  over  a  large  enough 
time  window.  Therefore,  the  program  also  has  the  option  that  if  ISUD=1,  then 
in  successive  10-minute  windows  the  noise  spectrum  is  recursively  updated 
frequency  by  frequency,  as  is  the  LTA  with  the  parameter  ISIGS  (analogous 
to  ISIG) ,  set  equal  to  -2  in  this  study.  As  we  shall  see,  with  this  modi¬ 
fication  the  optimum  filter  slightly  outperformed  the  bandpass  filter  when 
used  in  conjunction  with  the  Z-log  transform. 

An  option  exists  to  compute  the  "Z"  form  of  the  STA/LTA  detector.  If 
IZ=1  then  NDP  calculates  the  "Z"  variable  which  is  given  by  (LTA- STA) /vari¬ 
ance  where  the  variance  is  estimated  as  the  recursively  updated  value  of 
2 

(LTA-STA)  and  where  the  time  constant  is  taken  to  be  the  same  as  that  for  the 
LTA  itself.  If  IZL=1  then  the  log  of  the  STA  is  taken  before  the  LTA  and 
variance  are  computed. 

For  the  Z  detector  the  threshold  is  calculated  on  the  assumption  that 
the  Z  variable  is  normally  distributed  and  so  the  threshold  is  computed  via 
subroutine  QUANTF  which  gives  the  threshold  in  units  of  standard  deviations 
for  a  fixed  probability  of  false  alarm. 

For  the  other  detectors,  when  the  STA  is  made  up  of  squared  filtered 
values  (IASQ=1),  then  the  STA  is  assumed  to  be  distributed  as  chi-squared 
with  2BT  degrees  of  freedom  where  T  is  the  length  of  the  STA,  and  B  is  the 
bandwidth  of  the  filtered  noise.  B  is  calculated  from  the  spectrum  of  the 
filtered  noise  by  dividing  the  total  power,  (with  the  maximum  of  the  smoothed 
spectrum  normalized  to  1),  by  the  folding  frequency.  The  spectrum  may  be 
further  smoothed  by  IBWS,  typically  set  equal  to  15,  before  the  bandwidth  is 
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computed.  This  second  smoothing  makes  the  maximum  more  stable  and  therefore 
stabilizes  the  estimate  of  the  bandwidth.  If  a  bandpass  filter  is  used  then 
the  integral  for  the  bandwidth  is  taken  over  the  smoothed  noise  spectrum 
between  the  limits  of  the  filter,  (FLBW,  FHBW) . 

Once  the  2BT  degrees  of  freedom  have  been  calculated  then  the  threshold 
is  calculated  using  an  inverse  chi-squared  routine,  CHINV. 

If  a  calculated  threshold  is  not  giving  the  desired  false  alarm  rate, 
then  it  may  be  multiplied  by  the  factor  THDDL.  This  factor  is  generally  in 
the  range  0.9  to  1.1. 
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DISCUSSION  OF  RESULTS 


In  Table  I  and  Figure  1,  we  see  details  on  the  results  of  running 
various  of  the  program  options  discussed  in  the  preceding  section  on  the 
Pinedale  detection  tape  which  was  put  together  in  the  manner  described  in 
Appendices  I  and  III.  Run  1  is  the  "benchmark"  discussed  in  Appendix  III 
and  in  the  comments  field  in  Table  I  we  see  an  indication  of  the  parameters 
specifying  this  run.  (In  general  terms,  the  parameters  for  any  given  run 
are  the  same  as  those  for  the  run  immediately  above  it  in  the  Table  except 
for  those  which  are  noted.  Also,  unless  otherwise  specified,  parameters 
will  have  the  values  indicated  in  the  previous  section  as  "usual  in  this 
paper.")  Thus,  the  benchmark  run  number  1  squares  the  data  before  forming 
the  STA,  has  a  short-term  average  of  1.8  seconds,  and  requires  3/3  values  of 
STA/LTA  over  the  detection  threshold  which  is  determined  by  x2  with  2BT 
degrees  of  freedom  and  adjusted  by  the  factor  0.85. 

To  best  get  an  overview  of  the  relative  performance  of  the  detectors,  it 
is  useful  to  plot  them  on  a  graph  presenting  both  false  alarm  rate  and  proba¬ 
bility  of  detection  or  threshold.  The  usual  presentation  is  false  alarm 
rate  versus  probability  of  detection  for  a  fixed  signal  to  noise  ratio; 
this  is  termed  the  receiver  operating  characteristic  (ROC).  However,  in 
this  study,  we  choose  another  display,  that  of  false  alarm  rate  versus 
relative  magnitude  threshold.  The  false  alarm  rate  and  criteria  for  detec¬ 
tion  are  determined  as  discussed  in  Appendix  II,  page  2,  paragraph  4.  The 
relative  magnitude  threshold  is  determined  first  by  adding  up  the  detections. 

A  corrected  number  of  detections  is  then  determined  by  subtracting  an  esti¬ 
mate  of  the  number  of  detections  which  are  probably  false  alarms.  This  is 
determined  approximately  as  the  false  alarm  rate  times  the  total  signal 
windows  in  which  a  detection  is  not  expected  to  occur.  These  signal  windows 
are  taken  to  be  those  in  the  3rd  and  4th  S/N  level,  plus  1/2  of  those  in  the 
2nd  S/N  level.  Thus,  for  Pinedale,  with  a  total  of  124  signal  windows,  the 
expected  number  of  false  detections  for  case  1  is  (2.5/4)  x  124  x  (.5/60)  x 
(F/HR)  =  3.81.  The  (0.5)  comes  from  the  fact  that  the  signal  window 
is  0.5  min.  The  F/HR  =  116/(9.5  x  124/60)  =  5.9.  Thus,  the  corrected 
number  of  detections  is  42  -  4  =  38. 

To  convert  this  into  a  relative  magnitude  estimate  we  first  note  that  the 
average  difference  in  S/N  between  S/N  levels  is  0.3  magnitude  units  since  the 


signal  levels  are  decreased  by  a  factor  of  2  for  each  step. 

Thus,  if  detector  1  detected  each  signal  to  one  lower  level  than 
detector  2,  we  would  have  no  hesitation  in  saying  that  detector  1  was  0.3  m^ 
units  better.  Now,  there  are  about  30  signals  per  magnitude  level  so  that 
we  are  immediately  tempted  to  say  that  if  detector  1  detects  one  more  event 
than  detector  2,  then  it  is  .3/30  =  0.01  units  better. 

This  is,  in  fact,  the  interpretive  approach  we  shall  follow  and  we 
may  further  justify  it  as  follows. 

Even  if  all  of  a  set  of  signals  and  noise  have  exactly  the  same  long¬ 
term  stochastic  S/N,  the  fluctuations  of  the  noise  process  guarantee  that 
there  is  neither  perfect  detection  nor  perfect  non-detection  at  a  fixed  false 
alarm  rate.  In  fact,  Blandford  and  Wirth  (1973)  showed  that  one  could 
characterize  a  particular  automatic  power  detector  as  if  it  were  an  infin¬ 
itely  sharp  detector  as  a  function  of  S/N  but  that  the  S/N  ratio  varied  with 
a  standard  deviation  of  0.092  m^  units.  Thus,  even  if  all  "true"  S/N  values 
are  exactly  the  same,  an  improvement  in  the  detection  threshold  will  result 
in  a  few  more  events  being  detected.  For  example,  in  the  case  where  for 

detector  1  the  threshold  is  0.15  m,  below  the  fixed  (S/N)  we  would  expect  to 
.  15  b 

get  QC^Tjpj)  x  30  =  1.5  events.  If  the  threshold  were  improved  using  detec- 

05 

tor  2  by  0.1  m^,  then  we  would  expect  to  get  x  30  =  9.3  events  for  a 

difference  of  7.8  events,  almost  equal  to  the  10.0  we  would  get  if  the  S/N 

/x 

Z(x)dx  where  Z(x)  is 

the  unit  normal).  ° 

In  fact,  we  do  not  expect  that  the  S/N  values  for  all  signals  in  an 
S/N  level  are  the  same.  The  signal  spectra  change  in  shape  from  event  to 
event,  the  signal  lengths  change,  thus  changing  the  S/N  in  a  fixed  window 
length,  and  the  noise  spectra  change  as  a  function  of  time.  Thus,  it 
seems  reasonable  to  assume  that  within  any  given  S/N  level  there  is  a  rather 
broad  distribution  of  "equivalent"  S/N  values  that,  when  convolved  with  the 
detector  response  for  fixed  S/N  as  discussed  above,  results  in  an  additional 
event  detected  for  each  .01  m^  improvement  in  threshold. 
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TABLE  I 


This  table  gives  the  parameters  and  principal  results  of  the  detection 
runs  performed  by  the  authors.  The  successive  columns  give  the  run  number, 
the  identifying  name  of  the  run  of  program  DPP  which  appears  on  the  computer 
printout,  the  month  and  day  of  1980  on  which  the  run  was  made,  and  four 
columns  giving  the  number  of  false  alarms  and  detections  in  each  of  the  four 
S/N  bands.  The  seventh  column  is  the  sum  of  the  preceding  four. 

In  the  comments  column  ,  the  only  parameters  which  are  given  are  those 
which  change  from  the  preceding  run.  0  stands  for  the  optimum  filter,  1-2 
Hz  or  1-1.5  Hz  signify  the  appropriate  bandpass  filter,  SM  signifies  the 
optimum  smoothed  filter,  x2  threshold  signifies  that  the  threshold  was 
determined  separately  for  each  10-minute  window  using  x2  statistics,  Z  and 
Z-log  runs  have  the  threshold  determined  by  the  unit  normal.  D  is  an 
adjustment  factor  to  a  computed  or  fixed  threshold  to  modify  the  false 
alarm  rate.  S  is  the  STA  time  in  seconds,  and  1/3  or  3/3  signify  the 
Q/Q'  parameters.  If  Q/Q'  =  3/3,  then  the  effective  value  of  the  STA  is 
-  5/3  times  the  stated  value.  V  is  the  fixed  turn  on  threshold  for  STA/LTA 
detectors.  The  remaining  parameters  are  self-explanatory  or  are  detailed 
in  the  text. 
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TABLE  I 
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FALSE  ALARM/HOUR 


Pinedale  Data 


In  Appendix  V  are  presented  plots  of  the  signal  windows  together  with 
the  original  signals  before  they  were  buried  in  th^  noise.  In  Table  I  we  see 
the  results  from  the  different  runs,  and  in  Figure  1  each  run  is  represented 
as  a  point  in  a  plot  of  corrected  false  alarm  rate  per  hour  versus  relative 
magnitude  threshold.  Zero  magnitude  threshold  corresponds  to  a  corrected 
number  of  detections  equal  to  46.  Run  number  1  is  the  benchmark  run  discus¬ 
sed  in  Appendix  III.  We  see  in  Figure  1  that  since  it  has  a  low  relative 
magnitude  threshold  for  a  typical  false  alarm  rate  that  it  is  actually  one 
of  the  poorer  detection  systems  studied  in  this  report.  (Even  so,  however, 
it  is  only  about  0.18  m^  worse  than  the  best.)  The  points  connected  by  lines 
are  pairs  of  runs  which  differ  only  in  the  threshold  setting.  Thus,  they 
define  a  curve  such  that  as  the  magnitude  threshold  improves  the  false 
alarm  rate  increases.  For  example,  the  line  connecting  runs  2  and  3 
represents  a  threshold  change  from  V=1.85  to  v=1.7  for  a  3/3,  1-2  Hz  filter, 
absolute  value,  1.8  second  STA  detector.  We  see  that  the  lines  drawn  do  not 
cross  each  other  and  are  fairly  parallel  so  that  on  an  empirical  basis  it  is 
possible  to  say  that  if  point  A  is  to  the  left  of  a  line  with  such  a  slope 
drawn  through  point  B,  then  the  line  through  A  represents  a  better  detector. 

We  may  first  discuss  whether  a  window  of  3,  10,  or  30  seconds  is  best 
for  detection.  One  might  anticipate  that  10  seconds  would  be  best  because 
inspection  of  the  signal  length  estimates  in  Table  II  of  Appendix  III  shows 
that  they  are  fairly  uniformly  distributed  between  3  and  30  seconds.  That 
10  seconds  is  best  may  be  seen  by  comparison  of  runs  8,  10,  and  11  which  are 
identical  except  that  they  are  for  3,  10,  and  30  second  STA  windows  respec¬ 
tively. 

These  are  optimum  filter  runs  without  smoothing  and  with  threshold 
determined  by  automatic  x2  calculations.  We  see  that  the  10  second  window 
is  about  0.05  m^  units  superior  to  either  of  the  other  two.  In  other  calcu¬ 
lations  we  have  found  if  we  knew  which  signal  was  coming  in  advance,  then 
using  exactly  the  right  window  does  convey  an  advantage  of  about  0.1  magnitude 
units.  In  the  absence  of  such  information,  then  the  proper  thing  to  do  is  to 
pick  a  middle  value.  Comparison  of  runs  24,  25,  and  26  show  the  same  effect 
although  the  absolute  difference  is  only  on  the  order  of  0.03  m^  units  instead 
of  0.05.  The  difference  is  probably  within  the  range  of  statistical  fluctua¬ 
tions  for  these  calculations. 
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One  may  Chen  ask  which  of  several  filtering  procedures  is  optimum. 

For  this  we  may  compare  runs  10,  13,  14,  and  15  which  use,  respectively, 

the  optimum  filter  for  the  case  S/N=l  at  the  maximum  value  of  S/N;  the 

optimum  filter  for  the  case  S/N  <<  1;  the  prewhitening  filter  1/N, 

optimum  if  S/N  »  1,  and  otherwise  known  as  a  prediction  error  filter  if 

2  2  2 

causal;  and  finally,  the  Weiner  filter  S  / (S  +  N  )  for  the  case  S/N=l  at 
the  maximum  value  of  S/N. 

From  inspection  of  Figure  1  it  is  easy  to  see  that  the  prewhitening 
filter  is  about  0.25  inferior  to  the  other  three  filters.  Of  the  remain¬ 
ing  three  the  optimum  S/N=l  filter  is  about  0.02  m^  superior  to  the  other 
two.  It  was  a  bit  surprising  that  the  Weiner  filter  performed  so  well; 
perhaps  for  cases  where  there  is  better  S/N  at  high-frequencies  where  the 
signal  spectrum  is  not  flat,  then  the  "extra"  S  in  the  Weiner  numerator  will 
degrade  the  Weiner  filters'  performance. 

Along  the  line  of  the  optimum  filter  there  is  also  the  question  of  the 
optimum  signal  spectrum.  Again,  in  practice,  one  would  wish  to  find  the 
spectrum  of  the  expected  observations.  Typical  values  of  t*  are  in  the 
range  0.1  to  0.5  and  for  this  run  we  have  chosen  0.3.  To  comoare  to  another 
choice  of  t*=0.0  which  might  be  suitable  for  at  distances  up  to  10° 
we  may  compare  runs  11  and  12.  Here  we  see  that  for  this  data  base,  which 
consists  predominantly  of  teleseismic  P  waves, the  t*=0  spectrum  is  about 
0.06  m^  units  worse  than  the  t*=0.3. 

For  small  events  the  optimum  source  spectrum  for  P  waves  is  no  doubt 
a  flat  one.  However,  Figures  2  and  25  in  Sobel  and  von  Seggern  (1978)  show 
that  for  the  111^=5  level  events  which  have  been  scaled  down  in  this  study  a 
corner  frequency  of  1.0  Hz  is  typical.  Thus,  it  seems  reasonable  to  try 
detection,  assuming  a  corner  frequency  of  1  Hz.  In  this  case  we  may  compare 
run  25  which  has  the  standard  flat  spectrum  (corner  frequency  around  20  Hz) 
with  run  28  with  a  1  Hz  corner  frequency.  We  see  that  it  is  less  than  0.01 
m^  superior,  an  undetectable  amount.  In  any  event  we  would  not  use  a  1  Hz 
corner  frequency  in  practice  since  we  are  trying  to  detect  small  events,  not 
big  ones. 


As  to  whether  1/3  is  superior  to  3/3  in  detection,  it  would  appear  that 
it  is,  for  3-second  windows.  This  may  be  shown  by  a  comparison  of  runs  5 
and  7  which  are  identical  3-second  detectors  except  in  that  respect.  We  see 
from  Figure  1  that  the  1/3  is  approximately  0.07  m^  units  better  which  is 
probably  a  reflection  of  the  general  superiority  of  "incoherent"  detectors 
to  voting  detectors  as  discussed  by  Wirth  et  al  (1976).  In  the  1/3  detector, 
a  single  3-second  window  is  tested,  whereas  in  the  3/3  detector,  a  1.8  second 
window  is  slid  along  in  increments  of  0.6  seconds  and  3  consecutive  detections 
are  required,  a  kind  of  voting  procedure  over  a  total  of  3  seconds. 

We  may  also  ask  if  taking  the  absolute  value  of  the  filtered  data  for 
the  STA  is  worse  than  the  theoretically  optimum  procedure  of  squaring  the 
data.  Comparison  of  runs  2  and  5,  and  18  and  19  in  Figure  1,  shows  that 
for  the  average  of  the  two  cases,  squaring  is  superior  by  about  0.01  units, 
but  one  could  not  reject  the  hypothesis  that  there  was  no  significant  dif¬ 
ferences,  the  same  conclusion  reached  by  Vanderkulk  et  al  (1965)  and  Husebye 
(1972). 

One  may  see  that  the  threshold  is  not  very  sensitive  to  the  precise 
recursive  filter  bandpass  by  comparing  runs  4  and  6  which  have,  respectively, 
a  1-2  and  1-1.5  Hz  filter.  It  would  appear  that  the  1-1.5  filter  is  about 
0.01  m^  units  superior  to  the  1-2  Hz  filter,  a  difference  which  is  probably 
not  significant. 

We  may  also  evaluate  the  "Z"  detector  as  discussed  by  Shensa  (1976).  We 
shall  evaluate  only  the  broad-band  version  in  which  the  filtered  data  is  trans¬ 
formed  into  a  normal  distribution  using  a  mean  (LTA)  and  variance  recursively 

2 

computed  value  of  (STA-LTA)  .  In  the  present  program  the  constants  for  the 
recursive  calculation  of  the  variance  are  the  same  as  for  the  calculation  of 
the  LTA.  The  threshold  is  calculated  using  a  subroutine  for  the  inverse  of 
the  normal  distribution  (QUANTF).  The  Z  transformation  is  performed  if  the 
program  parameter  IZ=1.  If  IZL=1  then  the  log  of  the  STA  is  computed  before 
the  remaining  statistics  are  computed.  By  comparison  of  bandpass  filtered 
runs  20  and  23  with  21  and  24,  we  see  that  taking  the  log  does  not  seem  to 
improve  the  threshold.  The  run  that  corresponds  to  the  same  processing  but 
using  a  fixed  threshold  on  the  simple  STA/LTA  is  run  number  18  and  we  see  that 
it  is  equal  in  capability  to  the  Z  transform  result,  suggesting  that  the 
transformation  does  not  further  enhance  detection  capability  when  used  in 
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conjunction  with  a  fixed  filter.  However,  as  we  shall  see,  the  transforma¬ 
tion  does  keep  the  false  alarm  rate  constant  (as  do  the  x2  statistics), 
and  it  seems  to  enhance  the  capability  of  the  optimum  filter  as  we  see  below. 

Run  number  25  uses  the  optimum  filter  and  a  10-second  STA  with  the 
variable  threshold  computed  from  x2  statistics  with  2BT  degrees  of  freedom. 
When  compared  to  the  optimum  filter  using  the  log  Z  transform  we  see  the 
transform  is  about  0.07  m^  units  better.  Thus,  we  conclude  that  the  Z 
transform  neither  helps  nor  hurts  detection  capability  for  fixed  filters,  but 
together  with  optimum  filters  the  z-log  transformation  offers  a  significant 
advantage.  The  reason  for  this  is  unknown  at  present. 

To  examine  the  false  alarm  rate  question  in  more  detail  we  give  in  Tables 
II-l  and  II-2  the  false  alarms  occurring  for  the  Pinedale  tape  in  10  consecu¬ 
tive  2-hour  windows.  Since  detections  in  30  seconds  of  each  10-minute  window 
are  classed  as  detections  and  not  false  alarms,  the  numbers  given  are  actually 
the  number  of  false  alarm  detections  in  1-hour  and  54  minutes. 

Table  II-l  is  for  run  number  18,  the  best  recursively  filtered  trace; 
a  fixed  threshold  operating  on  a  10-second  squared  STA  out  of  a  1-2  Hz  fixed 
recursive  filter  with  detection  declared  on  a  single  threshold  crossing. 

Assuming  that  the  number  of  false  alarms  in  a  window  is  a  Poisson 
process  with  parameter  A  equal  to  the  average  number  of  detections  over  the 
whole  tape  one  cannot  reject  the  hypothesis  that  the  false  alarm 
rate  is  constant  with  90%  confidence  except  for  the  first  window.  For  this 
case  one  or  zero  detections  is  expected  to  occur  in  10  tries  with  <  0.01 
probability;  thus  we  might  conclude  that  the  false  alarm  rate  _is  varying  for 
this  case.  In  Table  II-2  we  see  the  results  for  run  31  the  log  Z  transform 
coupled  with  the  optimum  filter.  Here  we  see  no  evidence  of  a  non-constant 
false  alarm  rate,  as  would  be  expected  from  the  design  of  the  detector. 

NORSAR  Data 

Analysis  of  the  NORSAR  data,  discussed  in  Appendix  II,  yields  generally 
the  best  result;  this  result  is  only  slightly  better  (>  0.03  n^)  than  the  same 
procedure  using  a  3-second  window  (16),  but  much  better  (*  0.1  b^)  than  a 
30-second  window  (22). 


TABLE  II-l 


FIXED  FILTER  PWY  (RUN  18) 


S/N 

1/2 

1/4 

1/8 

1/16 

F/2  HR 

i; 

F 

D 

F 

D 

F 

D 

F 
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1 

1 

1 
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0 

0 

1 
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0 

2 

0 

1 

0 

1 
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0 

0 

0 

3 

0 

1 

0 

1 

0 
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0 

0 

1  prob  <  1  in 

4 

1 

1 

1 

0 

2 

0 

0 

0 

times  <  .01 

5 

1 

1 

0 

0 

3 

0 

0 

0 

6 

1 

1 

0 

0 
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1 

1 

0 

14 

7 

0 
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5 
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1 

1 

1 

0 

8 

0 

1 

1 

0 

1 

0 

1 

0 

9 

1 

1 

1 
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2 

0 

1 

1 

15 

10 

0 

1 

1 

0 
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0 

1 

0 

11 
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1 
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2 

0 

12 

0 

1 

1 

0 

0 

0 

3 

1 

9 

13 

1 

1 

1 

1 

1 

0 

1 

0 

14 

2 

1 

1 

0 

2 

0 

0 

0 

15 

0 

1 
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1 
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0 

2 

0 

16  prob  >  16  in 

16 

2 

1 
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1 

0 

1 

0 

0 

times  >  0.26 

17 

0 

1 

0 

1 

0 

0 

0 

0 

18 

1 

1 
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0 

0 

0 

2 

0 

5 

19 
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1 

0 

0 

0 

0 

0 

0 

20 

1 

1 

0 

1 

2 

0 

2 

0 

21 

2 

1 

1 

0 

0 

0 

1 

0 

9 

22 

0 

1 

0 

0 

0 

0 

0 

0 

23 

0 

1 

3 

1 

1 

1 

0 

0 

~  /. 
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1 

2 

0 

1 

0 

1 

0 

8 

25 

2 

1 

1 

0 

1 

0 

1 

0 

2b 

0 

1 

1 

1 

1 

0 

1 
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27 

0 

1 

0 

0 

0 

0 

0 

0 

8 

28 

2 

1 

0 

1 

0 

0 

0 

0 

29 

0 

1 

1 

0 

0 

0 

2 

1 

39 

1 

1 

0 

1 

2 

0 

1 

0 

9 

31 
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1 

1 

0 

0 

0 

0 

F/D  F/HR/D  ] 

c 

19 

31 

24 

14 

27 

5 

24 

3 

94/53  4.7/50 

i 

1 


-1 
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\ 

\ 
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t 
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TABLE  II-2 


OPTIMUM  FILTER,  LOG  Z  PWY  (RUN  31) 


S/N 

1/2 
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F/2  HR 

F 

D 

F 

D 

F 

D 

F 

D 

1 

0 

1 

0 

0 

0 

0 

1 

0 

0 

1 

0 

1 

2 

0 

1 

0 

3 

0 

1 

0 

1 

0 

0 

0 

0 
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4 

1 

1 

1 

1 

1 

0 

0 
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5 

0 

1 

0 

0 

1 

0 

0 

0 
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0 

0 

3 

0 

0 

0 

7 

7 

1 

1 

1 

0 

0 

0 

0 

0 

8 

1 

1 

1 

1 

1 

0 

0 

0 

9 

0 

1 

0 

0 

0 
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1 

1 

6 

10 

0 

1 

0 

0 

0 

0 

0 

1 

11 

0 

1 

1 

1 

1 

0 

1 

1 

12 

0 

1 

1 

0 

0 

0 

2 

0 

6 

13 

0 

1 

1 

1 

0 

1 

1 

0 

14 

1 

1 

0 

0 

0 

0 

1 

1 

15 

1 

1 

1 

1 

2 

0 

1 

0 

9 

16 

1 

1 

0 

1 

2 

1 

0 

0 

17 

0 

1 

0 

1 

0 

0 

0 

0 

18 

3 

1 

0 

0 

2 

0 

1 

1 

9 

19 

0 

1 

1 

0 

1 

0 

0 

1 

20 

1 

1 

0 

1 

3 

0 

0 

1 

21 

1 

1 

0 

0 

1 

0 

1 

0 

9 

22 

1 

1 

0 

1 

0 

0 

1 

0 

23 

0 

1 

0 

1 

1 

1 

1 

0 

24 

1 

1 

1 

0 

1 

0 

0 

0 

7 

25 

0 

1 

0 

1 

0 

0 

2 

0 

26 

1 

1 

0 

1 

1 

0 

r 

0 

27 

0 

1 

1 

0 

0 

0 

0 

0 

7 

28 

0 

1 

0 

1 

0 

0 

0 

0 

29 

1 

1 

1 

0 

1 

0 

2 

1 

30 

1 

1 

0 

1 

0 

0 

1 

0 

_ 1 _ _ _ 

31 

1 

1 

1 

0 

0 

0 

1 

0 

F/D  F/HR/DC 

17 

31 

12 

16 

24 

3 

20 

8 

73/58  3.7/56 
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Table  I  for  a  discussion  of  Table  format. 


RELATIVE  MAGNITUDE  THRESHOLD 


2.  False  alarms  versus  relative  magnitude  threshold  for  different 
detectors  on  NORSAR  test  tape. 


TABLE  IV- 1 


FIXED  FILTER  NORSAR  (RUN  24) 


S/N 

1/2 

1/4 

1/8 

1/16 

F 

D 

F 

D 

F 

D 
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D 

F/2  HR 
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0 

1 
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0 
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0 
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1 

0 
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.0 

1 

1 
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1 

0 
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.  9 
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1 
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0 

0 

9 
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7 
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0 

23 

0 
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1 
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1 
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0 

0 
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1 
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20 
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TABLE  IV- 2 


OPTIMUM  FILTER,  LOG  Z,  NORSAR  (RUN  18) 


1/2 

F  D 


1/4 
F  D 


0  0 
0  1 
0  0 
0  0 


0  0 


0  0 


1/8 

F  D 


0  0 


1/16 
F  D 


F/2  HR 


prob  9  or 
more  in  10 
times  > 

0.9 


prob  3  or 
less  in  10 


34  11 

0 

1 

1 

1 

U  1 

F/D 

F/HR/DC 

16  27 

17 

20 

14 

11 

21  5  68/63 

3.1/60 
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The  standard  "benchmark"  run  (1)  is  about  0.2  worse  than  the  optimum 
due  in  part  to  the  instability  of  the  optimum  filter. 

Comparing  runs  2  and  4  there  is  again  no  detectable  advantage  to  squaring 
over  taking  the  absolute  value,  despite  the  theoretical  result  that  squaring 
is  the  preferred  procedure. 

The  best  detector  which  does  not  use  the  Z  log  transform  is  run  21  which 
uses  the  fixed  1.5-4  Hz  filter,  a  10-second  windov ,  squaring  and  a  2BT  thres¬ 
hold.  It  is  about  0.05  m,  units  worse  than  the  best. 
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COMPARISON  WITH  WORK  BY  OTHER  CONTRACTORS 


In  Figures  1  and  2  we  have  presented  data  points  from  several  other 
detectors,  human  analysts  (B,  K) ,  the  MARS  process  (M),  Farrell  et  al  (1980), 
the  Texas  Instruments  power  detector  (T) ,  the  deflection  detector  (D) ,  and 
the  Walsh  detector  (W) ,  in  those  same  figures  the  Z-log  detector  operating 
on  a  10-second  window  as  discussed  in  the  preceeding  section  is  designated 
by  the  letter  Z.  Detections  by  these  processors  are  plotted  on  the  signal 
windows  in  Appendices  V  and  VI. 

In  this  section,  we  briefly  describe  each  of  these  detectors  and  make 
some  comparisons  between  them. 

The  results  for  the  human  analysts  (Richard  Baumstark  and  Raymond  Kimmel, 
with  over  20  years  AFTAC,  LASA  and  NEP  analysis  experience  between  them)  were 
obtained  as  follows.  The  data  were  plotted  out  at  develocorder  scale  on  a 
Calcomp  (the  same  scale  as  seen  in  Appendices  V  and  VI)  only  one  trace  on 
the  page,  in  windows  of  random  lengths  ranging  from  5  to  10  minutes.  The 
program  ensured  that  a  signal  did  not  begin  in  the  first  or  last  30  seconds 
of  a  sheet.  The  analysts  were  informed  of  the  method  of  plot  construction 
and  were  asked  to  analyze  the  plots  In  random  order  so  that  they  could  not 
"predict"  by  patching  plots  together  where  a  signal  would  occur.  The 
analysts  were  told  that  there  would  be  one  or  zero  true  signals  per  plot,  and 
that  they  should  try  to  make  fewer  than  three  calls  per  plot.  The  analysts, 
generally  speaking,  found  the  plots  difficult  to  analyze  but  thought  that 
the  presentation  was  "fair"i.e.,  not  biased  against  them,  except  possibly 
for  the  high-frequency  noise  at  PWY  which  comes  from  the  transformed 
quantization  noise.  In  Figure  1  the  letters  B  and  K  show  that  the  analysts 
liked  to  operate  at  a  very  high  false  alarm  rate  so  that  although  they 
detected  as  many  true  signals  as  tne  best  detectors,  their  high  false  alarm 
rate  leads  to  an  estimate  that  they  are  0.2  to  0.3  m^  units  worse  than  the 
best  detector  at  a  fixed  false  alarm  rate.  They  appear  to  be  close  in 
capability  to  the  prediction  error  filters,  a  result  in  agreement  with 
von  Seggern  (1977)  who  studied  prediction  error  filters  to  see  if  they  were 
suitable  as  first  motion  detectors.  After  the  analysts  had  analyzed  the 
NORSAR  plots,  they  were  asked  to  return  to  the  PWY  plots  and  reduce  their 
false  alarm  rate  by  discarding  those  detections  they  were  less  sure  of.  This 
resulted  in  substantial  improvement  in  performance  (see  B',  K'  in  Figure  1), 
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simultaneously  reducing  the  number  of  false  alarms  and  increasing  the  number 
of  "true"  detections.  (Of  course,  the  actual  number  of  detections  decreased 
by  (9,  10);  however,  this  was  overcome  by  the  reduction  in  estimated  false 
detections  due  to  the  decrease  in  the  false  alarm  rate.)  This  change  in 
performance  is  a  reverse  slope  to  the  trend  of  performance  of  every  other 
detector  discussed  in  this  paper  as  the  threshold  is  varied  and  suggests 
that  the  human  qualitatively  changes  the  type  of  processing  as  his  "thres¬ 
hold"  changes  and  does  not  simply  look  for  a  larger  value  of  some  parameter. 
For  example,  the  analyst  may  require  "duration"  in  addition  to  "frequency 
change"  if  he  wants  to  be  more  certain  instead  of  looking  simply  for  a 
greater  value  of  "frequency  change." 

In  Figure  2  for  NORSAR,  the  analysts  are  in  the  midst  of  the  performance 
points  for  other  detectors  and  in  fact  one  performance  is  as  good  as  the  very 
best  automatic  detector.  In  this  case,  the  analyst  reports  near-ideal 
peaceful  conditions  at  home  for  the  analysis,  and  that  he  was  able  on  this, 
and  some  other  occasions,  to  get  in  the  groove  of  the  analysis. 

These  results  suggest  that  under  routine  operational  conditions  with 
develocorder  film,  time  for  accommodation,  a  peaceful  environment,  and 
interest  in  the  work,  a  trained  human  analyst  can  be  at  least  the  equal  of 
an  optimum  detector  in  pure  detection  of  teleseisms.  However,  these  con¬ 
ditions  are  severe  ones  and  in  production  one  might  expect  the  machine  to 
outperform  the  analyst.  It  is  important  to  note,  however,  that  pure  detec¬ 
tion  is  only  a  part  of  the  analyst's  job.  They  must  also  discard  detections 
of  noise,  e.g.,  instrument  noise,  lightning,  electric  discharges,  glitches, 
machinery  noise,  locals,  etc;  and  must  detect  in  the  presence  of  other  events. 
Thus,  in  a  practical  sense,  operations  within  0.1  of  the  optimum 
as  indicated  in  Figures  1  and  2  for  the  "on-line"  system  is  generally 
satisfactory  until  the  serious  practical  problems  mentioned  have  been 
solved. 

This  result,  that  automatic  detectors  are  qualitatively  equal  in  pure 
detection  to  the  analyst,  has  been  the  general  understanding  of  the  SDAC 
analysts  who  worked  on  LASA,  and  on  NEP  after  the  DP  software  bugs  had  been 
fixed. 
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Examination  of  the  plots  in  Appendices  V  and  VI  show  that  the  analysts' 
time  picks  are  generally  more  accurate  than  those  of  the  other  detectors. 
(However,  of  the  plotted  detectors  only  the  MARS  processor  exercised  an 
option  to  compute  a  refined  start  time — which  was  plotted  in  the  Appendices— 
thus  strictly  speaking,  detection  time  comparisons  should  not  be  made.)  In 
fact,  the  analysts  so  often  picked  a  few  seconds  early  that  10  seconds  in 
front  of  the  signal  was  allowed  to  count  as  a  detection  on  the  theory  that 
the  analysts  detected  late  and  "backed  up"  to  an  incorrect  start. 
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The  TI  power  detector  (T)  is  based  on  results  provided  to  us  by 
J.  S.  Mott  of  ENSCO,  Florida.  The  detector  is  basically  that  of  a  Z-log 
detector  with  a  fixed  filter.  The  bandpass  filter  for  PWY  is  0.6  to  2.0  Hz, 
and  the  STA  length  is  9.6  seconds.  It  is,  therefore,  most  similar  perhaps 
to  our  PWY  run  21,  and  in  fact,  it  seems  to  lie  near  the  same  detection 
curve.  The  NORSAR  filter  is  a  0.8  Hz  high-pass  filter  followed  by  a  6  dB / 
octave  high  pass.  The  combination  is  down  by  6  dB  at  1.5  Hz  and  so  is 
comparable  to  our  1.5-4  Hz  filter.  It  is  similar  to  our  run  25  and  does 
indeed  seem  to  lie  on  the  same  detection  curve. 

The  maximum  deflection  detector  (D)  is  a  qualitatively  different  type 
of  detector  also  run  by  ENSCO  in  Florida.  In  this  detector,  after  highpass 
filtering  to  remove  the  microseisms,  an  FFT  is  taken  of  successive  3.2  second 
windows  and  the  window  is  slid  one  second  at  a  time.  Selected  individual 
periodogram  values  in  the  signal  "frequency  window"  (.93-3.86  Hz  for  PWY,  .93 
to  2.48  Hz  for  NORSAR)  are  tracked  as  a  function  of  time  by  incorporating  them 
into  the  "Z"  variable  form.  A  high  threshold  of  5.0  is  set  on  each  Z  value 
and  if  for  any  single  frequency  there  is  a  detection  then  a  system  detection 
is  declared. 

As  discussed  in  Appendix  I,  this  is  a  type  of  voting  detector  (Wirth 
et  al,  1976),  which  on  theoretical  grounds  is  expected,  all  other  things 
being  equal,  to  be  inferior  to  an  appropriate  broadband  detector.  Further¬ 
more,  by  allowing  a  detection  if  only  1  out  of  N  frequencies  detect  instead 
of  K  out  of  N  with  K  >  1  it  is  probably  not  even  an  optimum  voting  detector. 
(The  MARS  detector  is  similar  to  the  deflection  detector,  has  K=5,  and 
shows  better  performance.) 

In  fact,  we  see  in  Figures  1  and  2  that  the  deflection  detector  is  one 
of  the  least  successful  detectors,  although  it  does  out-perform  the  human 
analyst  in  some  cases.  This  points  up  the  difficulty  of  detailed  calibra¬ 
tion  of  automatic  detectors  by  comparing  them  to  analyst  picks  on  real  data. 

The  MARS  processor,  Farrell  et  al  (1980),  first  takes  in  a  window  of  about 
100  seconds;  they  are  tapered  by  10%  and  Fourier  transformed.  The  frequency 
range  of  0.25  to  5.0  Hz  is  split  into  20  bands  with  spacing  approximately 
equal  to  ~  5./20  =  .25  Hz.  The  width  of  the  filter  applied  to  each  band  is. 
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however,  ~  .08  Hz  so  chat  the  time  constant  of  such  a  filter  is  approximatley  1/.08 
=  12.5  seconds.  Thus  MARS  is  operating  with  a  -  10  sec  integration  time  as  was 
found  to  be  suitable  in  our  earlier  analyses.  The  envelope  of  the  narrow  band 
is  calculated  and  statistics,  mean  and  standard  deviation,  are  kept  of  the 
distribution  of  the  maxima.  If  5  bands  exceed  1.8  standard  deviations  in 
2.8  seconds,  a  detection  is  declared.  The  result  of  this  voting  type 

detector  are  seen  in  Figures  1  and  2.  For  PWY,  its  performance  was  sur¬ 
passed  only  by  the  Z  detector,  although  many  other  detectors  are  close  to  it. 
However,  for  NORSAR  it  is  surpassed  in  performance  by  many  others.  Perhaps 
this  is  due  to  the  greater  bandwidth  of  the  NORSAR  data  where  voting  detectors 
may  be  at  a  greater  disadvantage. 

Finally,  we  may  discuss  the  Walsh  detector.  In  this  case,  there  is  a 
digital  recursive  prefilter  to  the  data  cutting  out  the  microseism  energy 
below  1  Hz.  (This  is  required  because  of  the  large  sidelobes  of  the  Walsh 
transformation.)  Thereupon  the  Walsh  transform  is  made  of  3.2  second  data 
windows  and  the  Walsh  coefficients  prewhitened  with  the  average  from  the 
first  15  minutes  of  noise  and  summed.  This  value  is  compared  with  a  multiple 
of  a  selected  point  on  a  cumulated  histogram  and  detections  declared.  Thus, 
this  is  essentially  a  power  detector  following  a  high  pass  prefilter.  The 
algorithm  works  very  well.  The  requirement  of  recursive  frequency  prefilter, 
however,  means  that  its  computational  requirements  would  be  comparable  to 
those  of  other  standard  detectors,  all  other  system  and  programming  language 
considerations  being  equal. 
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SUMMARY 


As  discussed  in  the  abstract,  most  detectors  can  be  modified  to  perform 
in  a  near-optimum  fashion.  Probably  the  log-Z  plus  optimum  filter  detector 
and  its  future  enhancements  will  always  be  slightly  (~  0.05  m^)  better  than 
the  others  since  it  is  based  on  a  well-understood  statistical  foundation. 

The  subject  of  array  detection  is  equally  well-understood  and  has  been 
implemented  in  practice  at  LASA,  NORSAR  and  TFO.  Future  detection  research 
should  be  concentrated  on  experimental  studies  of  network  detection  and  on 
post-processors  which  determine  the  nature  of  the  received  signal.  For  this 
work  also  there  should  be  test  tapes  constructed  as  a  test  bed. 

The  comparability  of  the  automatic  detector’s  and  the  human  analyst's 
thresholds  points  out  how  hard  it  is  to  determine  "truth"  against  which 
automatic  detection  may  be  checked,  in  a  steady  data  stream.  Nonetheless, 
the  analyst  is  the  desired  judge  of  a  detector  in  an  absolute  sense.  It  is 
he  who  notices  in  the  course  of  event  analysis  that  the  detector  does  not 
miss  signals  he  would  have  called  and  it  is  he  who  notices  that  all  the 
detections  "show  something"  and  that  the  detector  does  not  deliver  up  to  him 
many  spikes  and  drop-outs.  To  achieve  this  status  is  mostly  a  question  of 
engineering  detail  and  not  of  theory.  By  these  criteria,  both  the  1971  LASA 
DP  and  the  current  SDAC  DP  are  successful  in  an  absolute  sense,  although 
relatively,  the  results  of  this  experiment  suggest  that  they  could  be  improved 
about  0.1-0.15  m,  . 
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INTRODUCTION 


In  most  detection  systems  of  which  I  am  aware  there  is  eventually  produced 
a  single  time  trace  on  which  a  detection  must  be  made,  either  by  an  analyst  or 
computer  program.  (one  exception  to  this  rule  is  FKCOMB  which  searches  w-K 
space  for  an  F-statistic  maximum).  Since  there  is  still  debate  within  the 
seismic  community  on  the  best  way  to  perform  this  basic  process,  I  propose  that 
we  first  attack  this  problem,  before  proceeding  to  3-component  and  array  pro¬ 
cessing  problems. 

In  general  we  do  not  know  what  the  wave-shape  of  a  signal  is  going  to  be 

so  that  match  filters  are  impractical,  but  we  usually  have  a  better  idea  about 

the  signal's  amplitude  spectrum  S  on  either  observational  or  theoretical 

grounds.  The  noise  spectrum  N  can,  of  course  be  easily  measured.  Freiberger 

(1963)  showed  if  both  N  and  S  are  Gaussian,  the  Neyman-Pearson  filter  with  the 

2  2  4 

highest  probability  of  detection  for  fixed  false  alarm  rate  is  S/(N(S  +N  )  ). 

Note  that  for  large  S  this  becomes  1/N,  a  pre-whitening  filter,  while  for  small 

S  it  becomes  (1/N)(S/N),  a  pre-whitening  filter  weighted  by  the  signal-to-noise 

2  2  2 

ratio.  Note  also  that  in  neither  limit  is  this  the  Weiner  filter  S  / (S  +N  ) 

whose  goal  is  not  detection  but  best  least-square  estimation  of  the  waveform. 

The  match  filter,  optimum  for  completely  known  signal,  is  the  small-S  limit 

2 

as  is  reasonable,  i.e.  S/N  where  in  this  case  S  is  complex. 

Workers  at  IBM,  Vanderkulk  et .  al  (1965)  were  fully  aware  of  this  theory 
and  designed  their  filters  using  the  small  S  limit.  They  showed  that  under 
certain  realistic  assumptions  a  band-pass  filter  would  be  a  good  enough  match 
for  (1/N) (S/N)  that  less  than  .05  m^  units  would  be  lost. 

The  Neyman-Pearson  detector  further  consists  of  squaring  and  integrating 

over  time  T  the  output  from  the  optimum  filter.  (The  output  (STA)  in  the 

2 

absence  of  signal  is  distributed  as  X  with  2BT  degrees  of  freedom,  where  B 
is  the  equivalent  bandwidth  through  the  filter.  In  the  presence  of  signal  it 
is  distributed  as  non-central  x  )•  IBM  (1965)  investigated  the  effects  which 
rectification  instead  of  squaring,  and  exponential  weighting  instead  of  in¬ 
tegration  would  have  on  the  detection  threshold  and  concluded  that  the  differ¬ 
ence  was  less  than  .05  m^.  Husebye  (1972)  came  to  the  same  conclusion.  This 
is  illustrative  of  a  general  theorem  -  "most  plausible  detectors  are  near- 
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optimum".  With  these  insights,  IBM  determined  the  detection  band-pass 
to  use  at  LASA,  and  for  many  years  the  detector  worked  to  the  satisfaction  of 
seismic  analysts.  In  my  own  mind  the  basic  theory  of  seismic  detection  on  a 
single  channel  has  been  closed  since  1963  although  "tuning"  of  the  filters  and 
integration  times  to  special  situations  still  seems  to  be  a  subject  for  research. 
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A  BRIEF  SURVEY  OF  SELECTED  RECENT  RESEARCH  IN  SINGLE-CHANNEL  DETECTION 

Gyoystdal  and  Husebye  (1972)  compared  prediction-error  filters  to  band¬ 
pass  filters  and  concluded  that  a  band-pass  filter  could  always  be  found  which 
would  out-perform  the  prediction-error  filter.  This  is  not  surprising  since 
we  saw  in  the  Introduction  that  a  1/N  filter  (which  is  what  prediction  error 
filters  attempt  to  achieve)  are  optimum  only  for  large  S/N,  whereas  (1/N) (S/N) 
is  optimum  for  small  (S/N).  Further  work  in  this  field  for  detection  purposes 
should  determine  an  evolutionary  filter  of  the  form  (1/N)  (S/N)  for  a  suitable- 
a-priori  S. 

LaCoss  (1972)  noted  that  the  false  alarm  rate  (FAR)  fluctuated  at  NORSAR 

and  developed  rules  for  varying  the  threshold  as  a  function  of  the  STA  variance 

to  hold  the  FAR  constant.  The  cause  of  this  variation  was  that  the  noise  field 

would  become  more  or  less  peaked,  resulting  in  a  different  bandwidth,  B,;  and 

thus  in  the  Neyman-Pearson  context  the  numbers  of  degrees  of  freedom  2BT  of  the 
2 

X  distribution  -  changed  and  the  threshold  had  to  change  in  order  to  keep  the 
FAR  constant.  LaCoss  did  not  cast  the  theory  in  this  context  but  defined  a 
stability  parameter  as  the  ratio  of  the  STA  mean  to  the  STA  variance.  This  can 
be  seen  to  correspond  to  a  measure  of  2BT,  and  so  is  monotonically  related  to 
the  FAR.  Bungham  and  Husebye  (1974)  report  that  continuously  varying  the  thresh¬ 
old  with  an  empirically  determined  function  of  the  stability  parameter  was 
successful  in  stabilizing  the  FAR  and  resulted  in  a  .025  m^  improvement  in 
detection  threshold.  Shensa  (1977)  and  other  workers  at  Texas  Instruments 
squared  the  STA,  which  is  the  optimum  process,  and  transformed  the  output  into 
a  normal  distribution  with  a  standard  deviation  which  varied  with  time.  (Note 
that  this  transformation  is  theoretically  impossible  if  the  original  time-series 
are  Gaussian).  The  attempt  however  seems  to  be  successful  in  a  practical  sense. 
Secoy  (1978)  evaluated  this  detector  on  array  beams  and  found  it  to  be  satis¬ 
factory  (there  were  however,  problems  with  spikes  and  dropouts).  This  de¬ 
tector  may  be  said  to  be  as  "standard"  as  the  IBM  detector,  and  some  normali¬ 
zation  of  this  type  seems  to  be  a  useful  idea.  Blandford  (1974)  found  that  an 
F-detector  operating  on  a  beam  for  one  month  also  performed  in  agreement  with 
theory  and  had  no  trouble  with  spikes,  dropouts,  local  events  or  noise  bursts. 

Wirth,  Blandford,  and  Shumway  (1976)  investigated  the  optimum  metnods  of 
combining  detections  from  several  stations.  They  found  that  "voting  detectors". 
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e.g.  declaring  an  event  when,  say,  3  out  of  5  Alaskan  stations  detect;  are 
slightly  worse  than  "incoherent"  detectors  in  which  the  STA  values  from  each 
station  are  added  together  and  detection  is  made  on  the  sum.  Ringdal,  Husebye, 
and  Dahle  (1972)  had  found  the  same  results  experl  mentally  comparing  voting 
and  incoherent  detection  at  NORSAR.  These  results  are  relevant  to  single 
channel  detection  because  one  may  consider  multiple  frequency  bands  as  multiple 
independent  stations.  Thus  detectors  which  break  up  the  spectrum  into  multiple 
frequencies,  detect  on  each  one,  and  then  declare  final  detections  on  the  basis 
of  some  rule  should  always  be  out-performed  by  simple  band-pass  detectors.  It 
might  reasonably  be  thought  that  exceptions  to  this  rule  will  occur  when  the 
signal  spectrum,  S,  is  very  poorly  known.  However,  the  work  of  Wirth  et.  al. 
(1976),  and  by  implication  the  work  at  NORSAR  where  signal  variability  is  very 
high  across  the  array,  show  that  in  this  case  the  incoherent  detector  has  truly 
a  substantial  superiority. 

Shensa  (1977)  investigated  three  different  single-channel  detectors.  The 
first  was  the  deflection  detector  where  successive  spectra  are  taken  of  the  time 
series  and  the  amplitude  at  each  frequency  is  tracked  over  time.  The  distri¬ 
bution  of  each  frequency's  amplitude  is  tabulated  and  by  use  of  logarithms 

transformed  into  the  unit  normal.  (Of  course  if  the  original  data  are  Gaussian, 

2 

each  amplitude  is  distributed  as  x  with  2  degrees  of  freedom  which  has  a  long 
tail  and  looks  much  like  a  log-normal  distribution).  Then  if  any  frequency 
exceeds  a  threshold  an  event  detection  is  declared.  This  amounts  to  a  voting 
detector  with  a  1/N  pre-filter.  The  second  detector  investigated  by  Shensa  is 
a  kind  of  incoherent  detector  called  the  average  deflection  power  detector 
where  the  unit  normals  are  averaged  over  a  signal  band.  Because  of  the  normali¬ 
zation,  this  might  amount  to  detection  with  a  1/N  per-f liter  followed  by  a 
bandpass  filter,  the  same  process  as  that  devised  by  IBM  (1965),  and  might  work 

rather  well.  However,  because  of  all  the  transformations  carried  out  on  the 
2 

multiple  x  distributions,  this  equivalence  cannot  be  confidently  asserted. 

The  third  detector  is  that  of  simply  summing  the  raw  power  in  the  signal  band 
and  then  converting  to  a  unit  normal.  This  detector  is  the  same  as  that  immedi¬ 
ately  above  without  the  noise  whitening  and  is  equivalent  to  the  one  evaluated 
by  Secoy  (1978). 

Savino  et.  al.  (1979)  have  discussed  the  MARS  detecuor  in  which  multiple 
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narrow  passband  filters  are  calculated  and  coincidences  in  detection  are  sought 
on  the  Hilbert  transform  envelope  of  each  output.  This  is  similar  to  the  de¬ 
flection  detector  of  Shensa  (1977)  and  has  the  same  weakness  that  detection  on 
some  number  of  traces  declares  an  event;  thus  we  have  basically  a  voting  detec 
tor.  Another  weakness  is  that  the  envelope  in  effect  integrates  over  a  time 
window  of  length  1/Af  where  Af  is  the  filter  width.  There  is  no  particular 
reason  for  this  to  be  equal  to  the  length  of  the  signal  T;  over  which  one 
should  integrate  for  optimum  detection. 

Allen  (1978)  developed  a  detector  with  a  view  to  making  it  simple  enough 

that  it  could  be  implemented  in  a  microprocessor.  His  detector  works  on  a 

linear  combination  of  the  squared  signal  and  the  squared  signal  derivative. 

Obviously,  it  may  be  impossible  to  work  the  frequency  response  of  this  process 

2  2  k 

into  the  form  S/[N(N  +S  )  ]  so  that  it  cannot  be  an  optimum  detector.  In  parti¬ 
cular,  if  there  is  a  large  noise  peak,  such  as  the  microseism  band,  a  prefilter 
will  be  needed  so  that  the  6  db/octave  due  to  the  derivative  function  can  have 
a  chance  to  be  used  to  advantage. 

Another  detector  developed  with  an  eye  to  microprocessors  is  that  due  to 
Herrin  which  uses  the  Walsh  transform.  Here  too,  a  narrow  noise  peak  can  lead 
to  problems  because  of  the  side-lobes  of  the  Walsh  function  response  to  a  pure 
sine— wave.  Each  of  these  last  two  detectors,  however,  comes  with  an  impressive 
array  of  ad-hoc  procedures  for  eliminating  non-statistical  false  alarms  such 
as  dropouts,  noise  bursts,  and  spikes. 
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DETECTORS  RECOMMENDED  FOR  STUDY 


°  The  IBM  detector  should  be  evaluated  in  its  routine  form  and  while 

operating  in  an  optimum  fashion.  Squaring  should  replace  rectification. 

Signal  spectra  should  be  estimated  for  explosions  and  earthquakes  of 

various  sizes  for  the  P  and  Lg  phases,  and  combined  with  typical  noise 

2  2k 

spectra  to  produce  optimum  S/  N(S  +N  and  (,1/N)  (S/N)  filters.  The 
integration  time  T  should  be  adjusted  separately  for  regional  and  tele- 
seismic  distances  and  the  improvements  achievable  over  the  standard  detec¬ 
tor  reported.  As  a  final  topic,  N  might  be  updated  at  regular  intervals. 
Rapid  updating  in  the  presence  of  large  events  would  yield  optimum  detec¬ 
tion  of  mixed  events,  and  of  Lg  in  the  presence  of  coda. 

°  The  Z  detector  and  possibly  deflection  detectors  could  be  evaluated  as 
examples  of  detectors  with  fixed  false  alarm  rates  and  as  examples  of 
detectors  which  use  the  "Z"  transform. 

°  The  MARS  Detector  could  be  evaluated  as  an  example  of  a  voting  detector. 

o  The  Allen  and  Walsh  detectors  can  be  evaluated,  with  prefilters,  to  see 
how  these  non-optimum  detectors  perform.  If  they  are  close  to  the  other 
detectors  in  performance  then  their  simplicity  may  recommend  them  in 
those  situations  where  an  analog  prefilter  is  available  or  not  needed. 

o  The  Allen  Detector  should  be  evaluated  for  installation  in  intelligent 
line  interfaces  (ILI's)  since  it  has  capability  for  rejecting  spikes. 
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EVALUATION 


I  suggest  that  the  performance  test  for  these  detectors  be  made  up  of  real 
signals  buried  in  real  noise.  The  noise  would,  however,  have  been  modified  by 
scrambling  the  phase  so  that  any  real  weak  signals  would  be  dispersed  through¬ 
out  an  FFT  window  and  would  be  very  unlikely  to  cause  false  alarms. 

The  technique  for  creating  signals  buried  in  noise  would  be  as  follows: 
Noise  Generation 

°  A  simple  detector  is  run  over  the  data  and  all  signals  are  detected  with 
S/N_>  some  fixed  threshold. 

°  With  these  detections  at  hand,  an  analyst  brings  a  first  4096  pt  time 
segment  up  on  the  tektronix  screen.  Assuming  that  he  sees  no  detection, 
and  that  the  detector  also  detected  none,  he  initiates  a  program  which 
performs  a  10%  cosine  taper  on  the  data,  Fourier  transforms  the  results, 
assigns  random  values  to  the  phase,  re-inverts  to  the  time  domain, 
applies  50%  cosine  taper  and  stores  the  result  in  core. 

°  Skipping  2048  points  into  the  raw  data,  the  same  process  is  repeated, 
and  the  result  is  added  to  the  last  2048  points  of  the  previous  data. 

(In  the  overlap  region  note  that  the  sum  of  the  two  cosine  tapers  is 
1-0).  The  previous  2048  points  are  read  out  to  the  subset  save  tape. 

°  The  process  is  continued,  except  that  when  there  is  a  detection  or  drop 
out  in  the  window,  the  previous  time  period  is  used.  However,  a  new 
set  of  random  phases  is  generated  so  that  the  noise  does  not  repeat 
identically.  Make  subset  individual  seismograms  8192  points  long  -  i.e. 
13  min  39.2  sec  long  for  N0RSAR  and  6  min  49.6  sec  for  PWY. 

°  The  final  result  is  a  tape  full  of  noise  whose  spectrum  changes  as  does 
the  real  noise,  but  which  has  no  signals  in  it. 

Signal  Additions 

°  Signals  will  be  added  every  13  min  39.2  sec*  exactly  on  the  second. 
Signals  will  be  selected  by  scanning  automatic  detections  with  high 
signal-to-noise  ratios  on  the  same  channel  under  consideration.  One 
minute  of  signal  preceeded  by  one  minute  of  noise,  the  dividing  line 
its  to  be  chosen  by  the  analyst,  will  be  spun  off  to  subset  tape. 

*  Changed  for  tapes  delivered,  see  Appendices  II  and  III. 
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Signals  are  to  be  fully  representative  of  locals,  regionals  and  teleseisms 
with  some  special  cases  for  shots. 

0  To  add  the  signals  to  the  noise,  the  analyst  will  skip  12  min  39.2  sec 
into  the  noise  tape  and  add  the  signal  at  such  a  level  e  that  it  could 
be  easily  detected  at  its  maximum  by  an  experienced  analyst.  Then  the 
same  signal  will  be  added  in  13  min  39.2  sec  later  at  <5  =  1/2  this  ampli¬ 
tude,  and  then  again  at  1/4  and  1/8.  To  make  for  a  smooth  transition 
of  the  noise,  the  signal  window  before  adding  will  be  tapered  by  a  25% 
cosine  taper;  and  the  noise  tape  will  be  multiplied  by  one  minus  the  pro¬ 
duct  of  e  6  and  the  cosine  taper,  Cl  -  £  6  cos)  so  that  the  mean  square 
noise  level  does  not  increase  through  the  signal  window. 

Other  Features  of  the  Test  Tape 

°  We  will  put  in  a  few  cases  of  large  signals,  spikes,  and  drop-outs 
followed  by  small  signals  so  that  those  workers  who  wish  to  test  their 
ability  to  detect  in  the  shadow  of  large  events  may  do  so,  and  so  that 
diagnostics  of  "glitch"  false  alarms  may  be  exercised  if  desired. 

°  We  suggest  24  hours  of  data  for  2  channels  this  will  result  in  2  X  104 
test  signals,  52  different  ones  at  four  amplitude  levels  each. 

°  Two  channels,  one  with  high  microseisms;  near  the  coast,  and  one  with 
low  microseisms  would  seem  to  make  a  good  test. 

°  Regional  events  will  be  represented  by  both  a  P-wave  and  an  Lg  window 
since  the  two  phases  have,  generally,  greatly  different  spectra. 

Suggested  False  Alarm  Rates 

Different  false  alarm  rates  are  suitable  for  different  requirements.  If 
one  is  considering  an  isolated  single  element,  then  FAR»l/hour  might  be  suitable. 
If  one  is  considering  a  voting  detector  operating  on,  say,  the  Alaskan  network, 
then  a  single  element  FAR  of  1/minute  might  not  be  too  high.  We  have  found 
that  for  operating  our  on-line  DP  in  conjunction  with  NEIS  stations,  that  our 
automatic  association  program  is  not  badly  confused  at  a  FAR  of  1/(5  minutes). 

I  suggest  for  this  test,  that  we  run  at  two  different  rates  1/(10  minutes  and 
1/(30  minutes).  For  this  experiment  we  need  a  reasonably  high  rate  so  that  we 
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can  be  statistically  certain  that  we  are  operating  at  the  advertised  rate.  Only 
if  each  worker  in  the  experiment  runs  at  the  agreed-upon  FAR  will  we  be  able 
to  compare  the  detectors'  performance.  Note  that  if  more  than  one  detector  is 
used,  e.g.  for  regionals  and  teleseisms,  then  each  detector  FAR  must  be  divided 
by  M  where  M  is  the  number  of  detectors. 

Suggested  Immediate  Time  Schedule 

By  February  6  ship  to  all  participants  4  hours  of  raw  data  from  the  two 
selected  channels  in  the  proper  format,  presumably  SDAC  subset.  Along  with 
this  data  we  will  send  a  list  of  analyst  picks.  These  data  should  enable  the 
various  participants  to  tune  up  their  programs  to  such  a  degree  that  when  the 
test  tape  arrives  it  can  be  run  through  immediately. 

By  March  15,  the  test  tape  as  discussed  above  can  be  ready;  and  within  a 
week  thereafter  we  could  produce  the  "blind"  tape. 

Suggested  Data  Channels 

For  the  data  channel  dominated  by  microseisms,  I  suggest  the  center  element 
of  the  C3  subarray  at  NORSAR.  We  have  several  years  of  continuous  data  from 
this  source  at  the  SDAC,  and  thus  there  are  plenty  of  explosion  signals  at 
regional  and  teleselsmic  distances  to  merge  into  the  noise.  The  data  is  very 
high  quality  and  is  conveniently  manipulated  with  software  in-house  on  both  the 
11/70  and  360/44.  The  only  drawbacks  to  this  data  are  that  it  is  at  lOsps  and 
that  there  is  a  sharp  cut-off  at  4.75  Hz.  This  is  a  slight  drawback  from  the 
point  of  regional  detection,  but  to  me  it  seems  that  the  many  advantages  out¬ 
weigh  the  drawbacks. 

For  the  data  channel  not  dominated  by  microseisms,  I  suggest  the  on-line 
data  which  we  receive  on  the  DPS  tape  from  Wyoming.  This  data  has  all  the  ad¬ 
vantages  of  the  C3  data  and  is  at  20sps.  I  know  of  no  drawbacks  to  this  data 
source;  one  other  advantage  that  it  has  over  the  C3  data  is  that  we  have  con¬ 
tinuous  3-component  data  also. 


A  few  comments  about  the  Alaskan  Data  in  comparison  to  the  NORSAR  data. 

In  general  the  Alaskan  data  does  not  have  the  high  dynamic  range  of  the  NORSAR 
data,  and  it  is  filled  with  spikes  from  line  discharges.  Furthermore,  it  is 
not  at  regional  distances  from  explosions  of  interest.  Thus,  despite  the  fact 
that  the  ALK  data  is  at  20sps,  it  seems  better  to  choose  the  NORSAR  data. 


SUGGESTED  FUTURE  WORK 


°  Follow  up  initial  tapes  with  tapes  with  signals  buried  at  random  as  a 
"blind"  test. 

°  Work  on  tapes  of  real  data  and  see  how  the  automatic  detectors  compare 
to  analyst  picks,  and  how  they  deal  with  spikes,  dropouts,  etc. 

°  Develop  tests  for  best  detectors  on  week-to-month  data  segments. 

°  Investigate  enhancements  possible  with  3-component  processing,  e.g. 

•  incoherent  detection  regarding  the  channels  as  independent 

•  detection  of  asymmetry  in  horizontal  components  with  Smart  detector. 
(The  F-detector  in  the  Smart  processor  is  not  reliable  due  to  con¬ 
tinuous  presence  of  signal,  the  processor  is  more  useful  to  deter¬ 
mine  azimuth). 

•  variations  on  Remode  (polarization  filters) 

°  Investigate  enhancements  possible  with  array  processing,  (a  comprehensive 
study  of  this  type  will  require  more  powerful  computers  or  an  array 
processor)  e.g. 

•  beamforming 

.  maximum-likelihood  detection 
.  F-detector 

°  Then  we  can  repeat  the  whole  show  with  LP. 
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APPENDIX  II 

JUNE  6,  1980  GEOTECH  MEMORANDUM  FROM  R.  R.  BLANDFORD, 
"TAPE  CONSTRUCTION  AND  BENCHMARK  DETECTION  RESULTS  FOR  NORSAR" 


MEMORANDUM 


WTELEDYNE  GEOTECH 

ALEXANDRIA  LABORATORIES 


TO: 

FROM: 

DATE: 

SUBJECT: 

i 
i 
i 

The  noise  data  are  from  24  May  1979  as  transmitted  by  satellite  at 
10  sps  from  the  center  element  of  the  C3  subarray  at  NORSAR  and  recorded 
at  the  SDAC.  This  day  was  selected  to  minimize  data  dropouts  and  system 
downtime.  The  data  were  plotted  as  a  check.  A  FORTRAN  computer  program 
was  written  to  create  "artificial"  noise  which  would  have  no  signals  in 
it.  The  program  operated  as  follows. 

A  4096  point  data  window  is  read  in  and  checked  for  signals  and  data 
dropouts.  If  there  are  none,  then  the  data  are  5%  cosine  tapered,  Fourier 
transformed,  the  phase  at  each  frequency  is  assigned  randomly,  and  the  data 
are  transformed  back  to  the  time  domain.  Then  the  noise  Is  tapered  with  a 
full  50%  taper  and  added  to  the  previous  window  (which  has  also  been  tapered) 
with  a  50%  offset.  Because  the  sum  of  the  two  50%  tapers  does  not  yield  a 
constant  root-mean-square  amplitude  this  fact  is  corrected  for  by  multi¬ 
plying  the  summed  time  series  by  an  analytical  function  of  time.  The 
successive  overlapped  portions  are  read  out  to  the  final  noise  tape  leading 
to  a  continuous  random  noise  field  whose  spectrum  and  amplitude  vary 
smoothly  as  does  the  true  noise  field. 

The  next  step  in  the  process  is  to  add  In  the  signals.  In  Table  1  we 

list  the  signals  used.  Signals  1  and  2  are  P  and  Lg  from  the  nuclear 

explosion  closest  to  NORSAR.  Figure  1  shows  the  P  signal  and  its  spectrum 
from  the  original  20  sps  high-rate  data.  This  is  the  raw  spectrum  uncor¬ 
rected  for  instrument  response.  Note  the  sharp  cut-off  at  5  Hz  such  that 
even  though  the  spectrum  is  flat  below  5  Hz  there  will  be  little  aliasing 
at  4  or  even  4.5  Hz. 

Seismograms  (3  and  4)  and  (5  and  6)  were  selected  from  Kazakh  test  sites. 
Signal  7  was  intended  to  be  from  Azgir,  but  due  to  a  mistake  no  signal  is  in 

the  data  window.  The  remaining  signals  are  selected  for  good  S/N  from  a  few 

day's  bulletins. 

The  signals  are  mostly  teleseismic  the  exceptions,  for  which  there  are 
good  signals,  being  signals  1,27,29,30,  and  31.  The  last  four  are,  unfor¬ 
tunately,  all  from  one  source  region,  Yugoslavia.  Because  the  other  signals 
are  teleseismic  and  since  they  had  to  have  large  S/N  in  order  to  be  suitably 
buried  in  the  noise,  they  will  have  lower  frequencies  than  regional  events 
and  generally  lower  frequencies  than  the  teleseismic  events  of  interest  for 
detection.  However,  corner  frequencies  scatter  greatly  as  a  function  of 
magnitude  so  that  although  the  sample  is  probably  biased  many  small  tele¬ 
seismic  events  would  have  these  spectral  shapes. 
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The  signal  windows  are  two  minutes  long  and  the  P  starts  were  centered 
in  the  window.  Thus,  there  is  1  minute  of  noise  in  front.  (Signals  2,  4, 
and  6,  being  Lg,  the  maximum  of  the  vertical  component  was  centered.)  The 
signal  windows  were  tapered  with  a  25%  cosine  taper  to  avoid  abrupt  starts 
and  stops,  then  added  to  the  noise.  Each  signal  was  added  to  the  noise  four 

times,  first  with  the  maximum  equal  to  1/2  the  raw  amplitude  maximum  in  the 

10  minute  window,  and  then,  at  1/4,  1/8,  and  1/16  the  amplitude.  Thus,  with 

the  34  signals  in  Table  I  there  are  a  total  of  4X34=136  10-minute  (6000  points) 

windows.  In  each  window  a  signal  begins  at  the  9th  minute  (point  5400).  The 
data  are  continuous  from  window  to  window. 

Due  to  an  error  which  occurred  during  the  signal  and  noise  merging  there 
is  a  graded  decrease  of  50%  in  the  rms  noise  in  the  30  seconds  preceding  the 
1/2  amplitude  signals.  With  an  LTA  of  -  30  seconds  or  less  this  makes  these 
signals  easier  to  detect  than  they  should  be.  However,  in  general,  the  first 
window  signals  are  easy  to  detect  by  analyst  or  macine  and  would  be  detected 
without  this  decrease  so  that  I  do  not  think  this  error  will  bias  the  use  of 
the  tape. 

To  provide  a  benchmark  we  have  run  a  standard  (optimum  frequency  filter, 
square,  LTA  ~  30  seconds,  STA  ~  1.5  seconds,  detect  in  3/3  successive  windows) 
detector  on  the  tape.  The  optimum  S/N^  filter  is  determined  from  the  first 
512  points  of  each  10  minute  window,  assuming  a  flat  source  spectrum  and  a 
t*  =  0.3.  The  resulting  filter  was  nearly  constant  from  window  to  window  and 
fairly  close  to  a  2-4  Hz  3rd  order  band-pass  filter.  Had  we  used  a  bandpass 
filter  then  the  detector  would  have  been  equivalent  to  the  standard  IBM 
detector  except  that  we  squared  instead  of  taking  absolute  values.  Processing 
was  10  times  real  time  on  the  360/44  -  most  of  the  time  was  taken  up  in  FFT's. 

A  detection  was  considered  valid  if  it  occurred  in  the  30  second  time 
interval  from  9  minutes  to  9  minutes  30  seconds.  At  the  false  alarm  rate 
(FAR)  of  one  every  (136x10) /161=8. 4  minutes  (see  Table  II),  we  expect  2.0 
false  alarms  in  each  S/N  gain  level  (2.0=(0.5x34)/8.4) .  Thus,  this  detector 
can  be  seen  to  be  operating  at  approximately  the  80%  or  better,  50%,  5%  and 
0%  detection  capability  at  the  four  gain  levels.  In  general  terms,  other 
experiments  using  absolute  values,  longer  STA  windows,  and  t*  values  ranging 
from  0  to  0.45  yielded  results  differing  only  slightly  in  detection. 

A  significant  improvement  would,  for  example,  be  detection  of  several 
signals  at  the  next  level  below  the  lowest  detected  level  for  each  signal. 

That  is,  detection  of  say,  5  of  the  underlined  signals  in  Table  II  without 
increasing  the  FAR  while  continuing  to  detect  the  higher  levels  would  be 
a  significant  improvement. 

In  general  terms  whenever  the  automatic  detector  triggered  the  signals 
could  be  detected  visually  by  me  when  the  raw  data  was  plotted  at  a  scale 
of  1"  to  5  seconds.  When  the  automatic  detector  could  not  detect  the  signal, 
neither  could  I.  It  might  be  that  a  more  extended  plot  scale  and  a  more 
experienced  analyst  might  detect  to  lower  levels;  however,  this  is  about  the 
characteristic  performance  of  the  SDAC  operational  11/70  detector. 
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Comments 


High  frequency,  5  sec 
Lg  1-4  Hz,  30  sec 

Clipped,  therefore  a  very  high  fre¬ 
quency  "event,  30  sec 
Lg  long  STA  would  work  well,  30  sec 

2. 5- 3. 5  Hz,  2  sec 

Lg  possibly  very  low  frequency 
No  signal  due  to  some  mistake 
2-4  Hz,  3  sec 
1-3  Hz,  5-10  sec 

1- 2  Hz,  2-10  sec 

2- 3  Hz,  2-5  sec 

1.5- 2. 5  Hz,  3  sec 
Nothing  visible,  2  sec 
1.5  Hz,  3  sec 

1  Hz,  3  sec 
1-2  Hz,  5  sec 

Very  short  signal,  2  sec 
1-2.5,  2  sec 

2  Hz  -  30  sec 

No  visible  signal,  3  sec 
1-3  Hz,  3-10  sec 
No  visible  signal,  too  weak 
1-4.5  Hz,  10-30  sec 

1  Hz ,  20  sec 

2  Hz,  2  sec 

3  Hz,  2  sec  with  precursor 
1-4  Hz,  20  sec 

No  visible  signal,  3-10  sec 
1-4  Hz,  20  sec 
1-4  Hz,  5  sec 

3- 4  Hz,  10  sec 

No  visible  signal,  nothing 
1-3  Hz,  5  sec 
1-5  Hz,  5  sec 


D  43/21-1  35/16  45/4  41/2 

Total  FA  =  164 

Expected  number  of  false  detections/gain  level 
-  802  -  50%  -5%  >0% 


2.0 


Comments  before  first  comma  are  obtained  by  inspection  of  the 
S/N=  1/2  traces,  comments  after  first  comma  are  obtained  by 
inspection  of  original  signals.  Comments  after  comma  were 
not  in  original  memo  of  6  June. 
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JULY  2,  1980  GEOTECH  MEMORANDUM  FROM  R.  R.  BLANDFORD, 

"TAPE  CONSTRUCTION  AND  BENCHMARK  DETECTION  RESULTS  FOR  PINEDALE" 
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MEMORANDUM 


WTELEDYNE  GEOTECH 

ALEXANDRIA  LABORATORIES 

TO:  Recipients  of  the  3  July  1980  VSC  Detection  Experiment  Tape 

FROM:  R.  R.  Blandford 

DATE:  2  July  1980 

SUBJECT:  Tape  Construction  and  Benchmark  Detection  Results 


The  noise  data  are  from  March  7  and  8,  1980,  as  transmitted  by  phone 
line  at  20  sps  from  Pinedale  Wyoming.  This  period  was  selected  to  minimize 
signals,  data  dropouts  and  system  down  time.  The  data  are  from  the  single 
vertical  KS  36000  instrument,  (channel  9).  A  FORTRAN  computer  program  was 
written  to  create  "artificial"  noise  which  would  have  no  signals  in  it. 

The  program  operated  in  the  identical  manner  as  that  discussed  in  my 
June  6,  1980  memorandum  about  the  30  May  VSC  10  sps  NORSAR  detection  tape. 

Table  I  gives  the  list  of  event  arrival  times,  distance,  latitude  and 
longitude  and  magnitudes.  It  can  be  seen  that  except  for  events  7,  8,  30  and 
31  the  data  are  teleseismic  and  of  large  magnitude.  In  general,  the  dominant 
periods  are  1  second  or  longer — much  longer  periods  than  at  NORSAR.  Signal 
envelope  lengths  range  from  3  to  50  seconds. 


The  signal  windows  are  two  minutes  long  and  the  P  starts  were  centered 
in  the  window.  Thus,  there  is  one  minute  of  noise  in  front.  The  signal 
windows  were  tapered  with  a  25%  cosine  taper  to  avoid  abrupt  starts  and  stops, 
then  added  to  the  noise.  Each  signal  was  added  to  the  noise  four  times,  first 
with  the  maximum  equal  to  1/2  the  raw  amplitude  maximum  in  the  10  minute  window, 
and  then,  at  1/4,  1/8  and  1/16  the  amplitude.  Thus,  with  the  31  signals  in 
Table  I,  there  are  a  total  of  4  x  31  ■  124  10-rainute. (12000  points)  windows. 

In  each  window  a  signal  begins  at  the  9th  minute  (point  10800) .  The  data  are 
continuous  from  window  to  window. 

To  provide  a  benchmark  we  have  run  a  standard  (decimate  to  10  sps,  optimum 
frequency  filter,  square,  LTA  30  seconds,  STA  1.8  seconds,  detect  in  3/3 
successive  windows)  detector  on  the  tape.  The  optimum  S/N2  filter  is  determined 
from  the  first  512  points  of  the  decimated  data  for  each  10  minute  window, 
assuming  a  flat  source  spectrum  and  a  t*  =  0.3.  Remaining  remarks  are  similar 
to  those  in  my  June  6  memo. 

In  Table  II  we  see  the  results. 


PWY  SIGNAL  LIST 


Sig. 

Order 

Seis . 
No . 

Sample 

Rate 

Date 

Arrival  Time 
—  1  min 

Lat . 

Long. 

A 

\ 

1 

i 

20.0 

12 

Sep 

78 

00:53:35.0 

29. 8N 

129. 5E 

89.6° 

5.3 

o 

2 

20.0 

13 

Sep 

78 

04:39:31.0 

26. 4N 

142. 0E 

84.8° 

5.0 

3 

3 

20.0 

15 

Sep 

78 

11:48:50.0 

48. 2N 

154. 3E 

63.2° 

6.1 

4 

4 

20.0 

16 

Sep 

78 

23:49:57.0 

25. 8S 

177. 9W 

92.8° 

5.1 

5 

5 

20.0 

25 

Sep 

78 

02:12:42.0 

41. 3N 

125. 5W 

11.9° 

4.8 

6 

6 

20.0 

25 

Sep 

78 

21:48:13.0 

50. 9N 

156. 2E 

60.7° 

4.6 

7 

7 

20.0 

24 

Jan 

79 

18:00:51.0 

37.  IN 

116. W 

7.5° 

4.5 

8 

8 

20.0 

23 

Mar 

79 

17:24:41.0 

26. 7N 

110 . 8W 

15° 

5.1 

9 

9 

20.0 

09 

Apr 

79 

02:03:40.0 

10. 5N 

82. 6W 

39.8° 

4.9 

10 

10 

20.0 

12 

Apr 

79 

00:25:03.0 

17. 8S 

178. 2W 

87.2° 

5.0 

11 

11 

20.0 

14 

Apr 

79 

02:51:34.0 

35. OS 

106. 8W 

77.5 

5.5 

12 

12 

20.0 

14 

Apr 

79 

10:11:27.0 

35. 7S 

102. 5W 

78.4° 

5.6 

(7 

sec  earlier 

13 

13 

20.0 

15 

Apr 

79 

06:31:15.0 

41. 9N 

19. IE 

83.8° 

6.0 

a  4.6  same 

loc) 

14 

14 

20.0 

15 

Apr 

79 

13:39:06.0 

3.4N 

82. 9W 

45.8° 

5.0 

15 

15 

20.0 

15 

Apr 

79 

14:54:34.0 

42.  IN 

18. 7E 

83.5° 

5.4 

16 

16 

20.0 

16 

Apr 

79 

10:16:09.0 

41. 9N 

19. 4E 

84.0° 

5.0 

1 7 

17 

20.0 

16 

Apr 

79 

21:03:46.0 

51. 2N 

179. 3W 

46.6° 

4.4 

18 

18 

20.0 

18 

Apr 

79 

13:27:49.0 

51. 4N 

170. 6W 

41.4° 

4.6 

19 

19 

20.0 

18 

Apr 

79 

18:09:44.0 

24. 4S 

67. 2W 

77.4° 

5.2 

20 

20 

20.0 

18 

Apr 

79 

22:15:38.0 

32. 4N 

41 . 2W 

53.8° 

5.1 

21 

21 

20.0 

21 

Apr 

79 

13:29:49.0 

52. 3N 

169. 6W 

40.6° 

4.5 

22 

22 

20.0 

22 

Apr 

79 

09:58:42.0 

32. 9N 

39. 6W 

54.6° 

5.4 

23 

23 

20.0 

26 

Apr 

79 

02:11:35.0 

33. 9S 

71. 9W 

83.7° 

5.4 

24 

24 

20.0 

28 

Apr 

79 

01:05:41.0 

18. 2S 

174.8W 

85.2° 

5.2 

25 

25 

20.0 

28 

Apr 

79 

11:49:18.0 

27. 5S 

71. 0W 

78.4° 

5.5 

26 

26 

20.0 

29 

Apr 

79 

14:07:35.0 

15. 2S 

178.5V) 

85.5° 

5.4 

21 

27 

20.0 

29 

Apr 

79 

16:10:24.0 

22. 5S 

177. 4W 

90.1° 

5.4 

28 

18 

20.0 

05 

May 

79 

20:12:38.0 

8. 4N 

71.  OW 

48.1° 

5.5 

29 

29 

20.0 

05 

May 

79 

20:16:18.0 

8.7N 

71. 1W 

47.8° 

5.3 

30 

30 

20.0 

06 

Aug 

79 

17:06:59.0 

37.  IN 

121. 5W 

10.5° 

5.4 

31 

31 

20.0 

15 

Oct 

79 

23:18:37.0 

32. 6N 

115. 3W 

11.0° 

5.7 
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1/2 


1/4 


1/8 


1/16 


Window 

Signal 

F 

D 

F 

D 

F 

D 

F 
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1 

4 

1 

2 

1 

1 

0 

2 

0 

1 

0 

5 

8 

2 

0 

1 

1 

1 

1 

0 

0 

0 

9 

12 

3 

0 

1 

1 

1 

0 

1 

0 

0 

13 

16 

4 

1 

1 

3 

0 

2 

0 

0 

0 

17 

20 

5 

2 

0 

1 

0 

1 

0 

0 

0 

21 

24 

6 

0 

1 

1 

1 

2 

0 

0 

0 

25 

28 

7 

0 

1 

2 

1 

0 

0 

2 

0 

29 

32 

8 

0 

1 

3 

0 

1 

0 

2 

0 

33 

36 

9 

0 

1 

0 

0 

2 

0 

1 

1 

37 

40 

10 

0 

1 

2 

1 

2 

0 

0 

1 

41 

44 

11 

0 

1 

0 

1 

0 

0 

2 

0 

45 

48 

12 

0 

1 

3 

0 

1 

0 

3 

0 

49 

52 

13 

0 

1 

0 

1 

3 

1 

1 

0 

53 

56 

14 

1 

0 

1 

0 

0 

0 

0 

0 

57 

60 

15 

1 

1 

0 

0 

1 

0 

0 

0 

61 

64 

16 

1 

1 

1 

l 

5 

0 

0 

0 

65 

68 

17 

0 

1 

0 

1 

2 

0 

0 

0 

69 

72 

18 

1 

1 

0 

0 

2 

0 

0 

0 

73 

76 

19 

0 

1 

0 

0 

1 

0 

1 

0 

77 

80 

20 

1 

1 

1 

0 

3 

0 

0 

0 

81 

84 

21 

3 

1 

0 

1 

1 

1 

0 

0 

85 

88 

22 

0 

1 

1 

1 

0 

0 

0 

0 

89 

92 

23 

0 

1 

0 

1 

1 

0 

2 

0 

93 

96 

24 

0 

1 

0 

0 

1 

0 

1 

0 

97 

100 

25 

1 

1 

2 

0 

2 

0 

1 

0 

101 

104 

26 

0 

0 

0 

0 

1 

0 

2 

0 

105 

108 

27 

0 

1 

0 

0 

2 

0 

1 

0 

109 

112 

28 

1 

1 

0 

0 

1 

0 

1 

0 

113 

116 

29 

1 

1 

1 

0 

0 

0 

1 

0 

117 

120 

30 

0 

1 

1 

1 

2 

0 

4 

0 

121 

124 

31 

0 

1 

1 

1 

1 

0 

3 

0 

17 

28 

27 

14 

43 

3 

29 

2 

Comments 

4  sec 
3  sec 

3  sec  clipped 
3  sec 

10  sec,  LF* 

3  sec 
30  sec  HF 
40  sec  LF 
3  sec 
3  sec 
30  sec 

20  sec,  v.  emergent 
30  sec,  clipped 
3  sec 
10  sec  LF 
10  sec 
3  sec  HF 
3  sec 

10  sec,  spikes  @  25 
3  sec  sec  in  noise 
3-20  sec 
20  sec 
40  sec 

5-10  sec  VLF 
10  sec,  10  sec  late 
25  sec 
3  sec 
12  sec 
3  sec 
~  60  sec 
-  sec 


*  LF  -  low  frequency,  HF  -  high  frequency 


TABLE  11.  False  alarms  (F)  and  detections  (D)  for  Pinedale  signals 
at  four  levels  of  S/N.  Detector  runs  on  decimated  data 
(10  sps) ,  optimum  S/N^  filter,  squared,  STA  *  1.8  sec, 
LTA  *  30  sec,  3/3  successive  detections.  False  alarm 
rate  1/10.7  min.  Expected  false  detections  per  S/N 
level  is  *  1.4. 
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APPENDIX  IV 

DESCRIPTION  OF  THE  SRC  (SDAC)  ON-LINE  DETECTION  SYSTEM 
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INTRODUCTION 


In  general  terms  the  on-line  detection  system  as  of  this  writing 
(January  1981)  represents  a  translation  to  FORTRAN  of  the  detection  algor¬ 
ithm  programmed  by  IBM  workers  in  microcode  and  machine  language.  For  the 
theory  of  this  detector  see  Vanderkulk  et  al  (1965),  and  for  the  earlier 
detailed  description  of  the  code  perhaps  the  most  accessible  reference  is 
the  SDAC  DP  Systems  reference,  several  copies  of  which  are  available  at  the  SDAC 
in  Alexandria.  For  all  except  the  most  determined  history  of  science  workers, 
however,  this  appendix  should  be  the  reference  of  choice  for  programming 
considerations.  The  present  on-line  program  has  features  not  available  in  the 
older  on-line  system  such  as  diagnosis  of  drop-outs,  complex  spikes,  and 
regionals.  The  regionals  are  diagnosed  using  the  spectral  ratio  proposed 
by  von  Seggern  (1977).  Another  new  feature  is  the  automatic  refinement  of 
start  times.  Because  the  original  IBM  system  also  beamformed  arrays  it  had  a 
feature  not  in  the  present  system,  i.e.,  after  detections  the  array  of  beams 
would  be  scanned  in  space  and  time  to  find  the  beam  with  the  largest  STA/LTA 
and  for  which  it  or  its  neighbor  had  been  the  largest  k  consecutive  times. 

That  beam  would  be  the  detection  beam. 

When  the  IBM  system  was  reprogrammed  in  machine  language  to  handle 
multiple  arrays  a  mistake  was  apparently  made,  resulting  in  inconsistent 
detection  performance.  To  remedy  this  problem,  the  system  was  reprogrammed 
in  FORTRAN  to  run  under  UNIX  on  the  PDP  11/70  at  the  SDAC.  To  provide  as 
much  continuity  as  possible,  the  program  was  made  functionally  as  identical 
to  the  IBM  system  as  possible. 

Because  of  the  use  of  microcode  it  was  possible  to  use  only  one  filter 

in  the  IBM  system.  Our  on-line  system  uses  a  filter  designed  according  to 
2 

S/N  principles  for  station  BFAK.  It  is  a  1.2-3  Hz  bandpass.  As  we  have 
seen  elsewhere  in  this  report,  this  is  not  optimum  for  the  on-line  PWY  channel. 

As  a  result  of  the  work  discussed  elsewhere  in  this  report,  the  on-line 
system  will  soon  be  changed  to  allow  different  filters,  to  square  the  data, 
to  detect  on  only  a  single  threshold  crossing,  and  perhaps  to  use  the  Z-log 
transform  and  perhaps  to  use  a  10-second  window  for  detection.  Other  work  is 
also  underway  on  incoherent  detection. 
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QUALITATIVE  DISCUSSION  OF  THE  DETECTION  SUBROUTINE  (STALTA) 

Subroutine  STALTA  contains  the  detection  algorithm.  Inputs  to  STALTA  are: 

CHANID  -  a  14-character  string  identifying  the  current  channel 

RAW  -  a  floating  point  array  containing  the  demultiplexed  and 

degained  data  for  the  current  channel  and  time  window. 

ICHAN  -  an  integer  array  containing  parameters  for  the  current 

channel 

CHAN  -  a  real  array  equivalenced  to  ICHAN 

XX  -  real  array  storing  previously  filtered  data  used  to 

eliminate  transient  ringing  in  the  recursive  filter 
STA  -  integer  array  containing  current  set  of  short  term  averages 

Y  -  integer  array  containing  intermediate  averages 

LTA  -  integer  value  of  current  Long  Term  Average 

QR  -  integer  array  containing  "detections"  for  'Q  out  of  Q' 

criterion 

DETN  -  logical  flag  indicating  "detection  in  progress" 

DISP  -  logical  flag  indicating  waveform  output  should  be 

produced 

DSTA.DLTA  -  floating  point  representations  of  "instantaneous"  STA  and  LTA 

DTIME  -  floating  point  start  time  of  detection 

STATUS  -  character  description  of  detection  -  'drop-out’, 

'local',  etc. 


1)  Given  input  data  array  of  30  seconds  worth  of  data  structured  as: 


7.5 

15 

H 

i - 

(1.2-3  Hz,  6-pole  Butterworth) 


recursively  filter  segments  2  and  3  onto  a  new  array  called  FILTER. 
Segment  3  is  saved  to  be  used  as  segment  1  of  the  next  pass.  Thus 

pass  (n-1) 


*  filtered 

|_J^I 

Pass  n  _ _ _ 

f i leered  XX 

pass  (n+1) 

filtered 

2)  Turn  on  and  turn  off  threshold  are  defined. 

3)  Begin  moving  average  over  points  taken  "S"  at  a  time.  This  corresponds 
to  sliding  STA 


trrn,t 


The  LTA  is  updated  whenever  the  STA  window  has  moved  into  a  completely 
new  set  of  points.  Thus,  after  STA  4  was  computed,  the  LTA  would  be  updated 
using  the  data  from  STA. 
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The  absolute  value  of  the  data  is  used  i.e., 

STAl  »  2\-f.l4crCO 
L*  t  1 

This  is  for  continuity  with  earlier  detector,  squaring  is  optimum. 


4)  The  STA/LTA  threshold  is  compared  against  the  current  threshold.  If  it 
exceeds  that  threshold,  a  flag  is  set. 


5)  If  STA/LTA  exceeds  the  threshold  for  Q/Q'  time  segments  (for  on-line  DP 
this  amounts  to  3  successive  exceedings) ,  we  may  declare  a  detection. 


6) 


£ 


See  decision  matrix  for  full  definition  of  start.  Opti¬ 
mum  Q/Q'  is  1/X,  any  x.  3/3  is  used  for  continuity. 


3 


Given  we  have  a  detection  start: 


a)  Compute  the  maximum  STA  of  the  three  used  in  declaring  the  event 

save  this  to  be  output  as  the  amplitude. 

b)  Refine  the  start  time  by  lowering  the  threshold  and  back  up 


t4o 


c)  Check  for  dropout 


1  ,3S<  3 

5a* p It  ro4e  n  to 
O^rriA+lwA 

r 

r  t-sur  K  -1 

we  check  the  RAW  (not  filtered)  data  for  ISR/2  consecutive 
zeros.  If  these  are  found,  we  assume  a  data  gap  and  bypass 
further  processing.  A  search  window  beginning  earlier  is 
desirable  to  check  for  data  re-acquisition.  A  large  number  of 
consecutive  zeros  are  needed  for  under-quantized  data  as  is 
characteristic  of  ALK.  For  NORSAR  two  zeros  in  a  row  would 
be  diagnostic. 


d)  Check  for  type  1  spikes  ^  see  spkchk 
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e)  FFT  computation  =  mean  removal,  10%  split  cosine  bell  taper. 

Inverse  X  form  after  removing  low  frequency  components 

f)  Check  for  spikes  on  filtered  data 

g)  Compute  periodogram.  Take  spectral  ratio  of  low  frequency  to  high 

frequency.  Categorize  as  "local"  or  "teleseismic"  via  spectral 
ratio 

End  of  Detection  Start  processing. 

If  Detection  End  -  i.e.,3  consecutive  STA's  less  than  end  threshold.  We  report 
everything  regarding  this  detection. 

LTA  is  updated  every  R  STA's. 

LTA  update  is  skipped  if  the  detection  was  a  dropout  or  spike. 

Two  different  LTA  decay  constants  are  used  depending  on  whether  a 
detection  is  in  progress  or  not.  This  is  to  close  down  a  detection  rapidly 
so  that  later  phases  will  be  detected. 

The  start  and  end  threshold  are  updated  to  reflect  the  new  LTA  . 


Recursive  Filter  Algorithm 
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P-1 


F(nAt)  =  £  ap*f ((n-p)At)  +  p£^  bp*F((n-p)At)) 

Letting  At  =  unity  yields  the  formally  simpler  expression 


P-1  P-1 

F(n)  *  I  a  *f (n-i)  +  £  b  -F(n-i) 

i=0  1  i=l  1 

where  f  represents  the  raw  data 

F  reprsents  the  filtered  data 
P  represents  the  degree  of  filter. 

Detection  Algorithm 

A)  Intermediate  STA  in  our  case  S  *  3  for  all  channels 

S-l 

y(n)  =  £  | F (n  -  c  | 

i=0 


B)  Short  Term  Average 


STA(n)  =  STA(n-l)  -  y  (n-R)  +  y(n) 

R  =  3 

c)  Long  Term  Average 

LTA(n)  =  e-°STA(n)  +  (1  -  e“°)  LTA(n-l) 
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STA 

>TH 

DETN 

Q/Q' 

No 

Arrival 

Start 

Continue 

End 

N 

N 

N 

X 

N 

N 

Y 

- 

- 

- 

- 

N 

Y 

N 

X 

N 

Y 

Y 

X 

Y 

N 

N 

X 

Y 

N 

Y 

X 

Y 

Y 

N 

- 

- 

- 

- 

Y 

Y 

Y 

X 

:  i 

Frequency  Domain  Spike  Check 


In  the  current  detector 
Q/Q'  is  3/3. 


Not  to  be  confused  with  deglitch  algorithm  which  is  designed  to 
eliminate  "single  point"  glitches. 


Deglitch  ->■  if  | f (i  -  1) |  <  G,  and  jf(i)|  >  10G,  and  |f(i  +  1) |  <  G 
f(i)  f(i  -  1) 


Normally  G  =  20.0 


Despike:  A  specific  algorithm  for  complex  ALK  spikes. 


Take  window  6.4  sec  -  same  as  FFT  window  of  64  points  at  10  Hz  of  data 
around  possible  detection 


a)  IF  i 


max 


-  i 


min  1 


<  10  AND 


(  |2f£|  >  1.5  OR  pS*|  <  i/i.s 
l  'min’  'min' 


THEN  define  detection  as  TYPE-1  spike 
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£  F(u>) 

b)  Apply  FFT  to  data;  compute  spectral  ratio  R  =  mf»2  Hz 

2  Hz 

]£  F(u) 
w=  1  Hz 


c)  Zero  transform  components  <  0.3  Hz  and 

inverse  transform  to  remove  microseism  trends  so  that  small  ALK 
spikes  may  be  analyzed. 

d)  for  new  data  window  apply  test  of  a)  above  but  with  1.5  replaced  by 
2  to  check  for  Type  2  spikes. 

e)  if  not  a  spike 

log  R  <  .2  "teleseism" 

log  R  >  . 2  "regional" 

if  "spike"  but  log  R  <  -.15  "teleseism". 


/ 

teleseisms  OK 

V 

V 

/ 

spikes  OK 

> 

locals  OK  \ 

\ 

must  be 
teleseism 

* 

-  0.15  log  R  +  .2 


Summary  Graph 


SELECTED  SUBROUTINE  LISTINGS 


RIT 


Feb  5  16:37  1981  rit.h  Page  1 


c 

c...  Detection  Criterion  Common 
c 

c...  NPTS  -  No.  of  DATA  points  to  process 

c...  S  -  No.  of  DATA  points  to  average  for  intermediate  sta  (Y) 

c...  NS  -  No.  of  Y's  to  average  for  STA 

c...  R  -  Interval  in  STA  units  between  LTA  update 

c...  Q  -  q  in  q/q'  criterion 

c...  QP  -  q'  in  q/q'  criterion 

c...  TRNON  -  Turn-on  Threshold 

c...  TRNOFF  -  Turn-off  Threshold 

c...  TSUBC  -  Normal  LTA  Decay  Time  constant 

c...  TSUBCD  -  Decay  Time  constant  during  detection 

c 

c...  IDX  -  Coherency  distance  for  Spkchk  routine 
c...  BRATIO  -  'Big  Spike'  ratio  "  "  " 

c...  S RATIO  -  'Small  Spike'  ratio  "  "  " 

c 

c...  I FLO  -  Low  FFT__channel  for  Spectral  Ratio  Calculation 
c...  IFMID  -  Middle  “  "  "  "  "  " 

c...  I FBI  -  High  "  "  "  "  " 

c...  IFCUT  -  Cutoff  FFT^_channel  for  inverse  filtering 
c...  SPECT  -  Local/Teleseism  Discrimination  Ratio 
c...  BCOEF  -  B-coef f icients  for  recursive  filters 
c 

Real  TRNON,  TRNOFF,  TSUBC,  TSUBCD 
Real  BRATIO,  SRATIO,  SPECT,  BCOEF(0:6,3) 

Integer*2  IDX,  NPTS,  R,  Q 

Integer*2  IFLO, IFMID, IFHI,  IFCUT 

Common  /rit/  TRNON,  TRNOFF,  TSUBC,  TSUBCD, 

&  BRATIO,  SRATIO,  SPECT,  BCOEF, 

&  NPTS,  R,  Q,  IDX,  IFLO, IFMID, IFHI,  IFCUT 

c 

c  End  of  rit.h 
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Dec  10  10:17  1980  sites. h  Page  1 


c 

c  •  •  < 

c 

c  •  •  < 

C  •  •  i 

c  •  • 

C  •  «  i 

c 


c 

c .  *  • 
c «  .  • 
c .  .  . 
c. .  . 
c. .  . 
c. .  . 
c. .  . 
c 


Channel  Parameters 
Parameter  (Maxcha“7) 

NDET  -  Number  of  channels  to  be  proccessed 
LTA  -  Current  LTA  value 
STA  -  ST A  array;  OP  points 

OR  -  Work  area  for  q/q'  calculation;  QP+1  points 

Integer*2  NOET 

Real  LTA(Maxcha),  ST A( 5, Maxcha) 

Integer*2  QB(5, Maxcha) 


XX 

DETN 

DLTA 

DSTA 

□TIKE 

STATUS 

CHANID 


Filter  save  area;  75  points 

Flag  indicating  'detection  in  progress' 

LTA  at  start  of  detection 
STA  at  "  "  " 

Time  at  "  "  " 

Alphanumeric  detection  type  (Local.  Spike^l,  etc.) 
Alphanumeric  channel  id 


c 

C  a  •  a 
C  •  •  • 
C  •  a  • 

c 


Real  XX(75, Maxcha) 

Logical  DETN(Maxcha) 

Real  0ETECT(6 .Maxcha) 
Character*10  STATUS(Maxcha) 
Character*14  CHANID(Haxcha) 


-  Site  number 

-  Displacement  to  channel  in  DPS  record 

"  "  "  in  site 

-  Increment  to  next  point  in  site 

-  Calibration  in  nm/count  for  RAW  data 

-  Delay  of  site  data  relative  to  DPS  Tape  TOD 

-  Seconds  of  Data  Processed 

-  Number  of  Local  Detections 

”  "  Teleseismic  " 

-  "  "  Spike- 1 

-  "  "  Spiked 

"  "  gram  files  produced 

Integer*4  I CHAN ( 12, Maxcha) 

Real*4  CHAN ( 12, Maxcha) 

Equivalence  (ICHAN<1,1),  CHAN(l,l)) 

Storage  for  Demultiplexed  data 

Note  that  the  last  150  elements  of  RAW1  are  not  used  in  the 
demux  operation  and  are  free  for  other  uses. 

Real  RAW0C150) 

Real  RAWK300, Maxcha) 

Real  RAW2( 300, Maxcha) 

Equivalence  (RAWlU  ,1)  ,RAW2(151 ,1)  ) 

Equivalence  (RAW0(1),  CHAN(l.l)) 
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c  •  •  • 

ICHAN  contents: 

c 

c  •  •  • 

1) 

Integer 

C  •  •  e 

2) 

Integer 

C  •  •  a 

3) 

Integer 

Cm  .  . 

4) 

Integer 

Case 

5) 

Real 

Cam* 

6) 

Real 

Ca  a  a 

7) 

Integer 

Can 

8) 

Integer 

Ca  a  a 

9) 

Integer 

C  a  •  • 

10) 

Integer 

C  a  a  a 

11) 

Integer 

C  a  a  a 

12) 

Integer 

i 


Dec  10  10:17  1980  sites. h  Page  2 


Common  /sites/  RAW2.RAW0,  XX, DETECT, 

&  NDET,LTA,STA,QR ,  CHANID,  STATUS,  DETN 

c 

c  End  of  sites. h 


* 
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Feb  6  13:50  1981  stalta.f  Page  1 

Subroutine  Stalta  (CHANID,  RAW,  ICHAN,  CHAN,  XX.STA, 

&  LTA,  OR,  DETN, DETECT, STATUS) 

Parameter  (ISR  “  10) 

Parameter  (S  ■  5,  NS  "3,  OP  "3) 

c...  STALTA  computes  the  short  &  long  term  averages  of  the 
c...  filtered  input.  The  ratio  of  STA/LTA  is  compared  to 
c...  a  threshold  value  to  define  a  detection, 
c 

Include  'rit.h' 

Include  'run.h' 

Character*14  CHANID 
Integer*4  ICHAN(*) 

Real*4  CHAN(*) 

Real  RAW(*) ,  XX(*),  STA<*).  LTA 

Integer*2  0W(*) 

Logical  DETN 
Real  DETECT(6) 

Character  STATUS*10 
c 

Real  calib,  delay 
Complex  xcl(64) ,  xc2(64) 

Real  dt,  df,  dtlta,  dtsta,  fltlta 
Real  qsum,  sumlo,  sumhi,  fmax 
Real  e2sig,  sigf,  e2sigd,  sigdf 
Real  Btatim,  rawtim.f iltim 
Real  xrl(64),  xr2(64) 

Real  filteq(300),  filter(150) 

Equivalence  (f ilter(l) ,f ilteq(76)) 

Real  tsrt.tend.th,  cursta 
Integer*2  i,j,  il,i2 
Integer*2  iknt,  nsta 
Integer*2  detptr.f ilptr.rawptr.f ftptr 
Logical  big,  qqp,  first 
Logical  detsrt,  detend,  saqout 
Logical  drop,  BSPK,  SSPK, valid 
c 

Data  first/ .true./ 

Data  filteq  1200*0.0/ 
c 

c...  Initial  Setup  Section 
c...  Done  only  if  LTA  “  -999.9 
c 

If  (LTA.eq.-999.9)  Then 
c 

c...  Define  scaling  factors 

C  dt  -  1.0/Float<ISfl) 

df  -  1 .0/ (64*dt) 
dtsta  ■  S*dt 
dtlta  -  (R*S)*dt 
e2sig  -  1.0  -  Exp(-dtlta/ tsubc) 
sigf  -  1.0  -  e2sig 
e2sigd  ■  1.0  -  Exp(-dtlta/tsubcd) 
sigdf  -  1.0  -  e2sigd 
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c 

c. .  . 
c 


20 


c 

c  •  •  • 
c 

c 

c  •  •  • 

C  •  •  • 

c 

40 

42 

c 

c. .  . 
c. .  . 
c 


c 

c  •  •  « 
c 


c 

c. .  . 
c. .  . 
c 


Compute  initial  LTA  and  thresholds 

Call  Rfil  (NPTS,BCOEF(0,ICHAN(12)),RAW(151), filter) 
LTA  -  0 

Do  20  i  -  1 ,S*NS 

LTA  ■  LTA+Abs(filter(i)) 

If  (first)  Then 
first  “  .false* 

Call  Header 
End  if 

saqout  *  .false, 
statim  “  timei+dtsta 
fltlta  -  Float( LTA) / (S*NS) 

DETECT (4)  -  0.0 
DETECT(5)  -  0.0 

Call  Report  (CHANID, statim, 'Initial  ' ,0 .0 , fltlta, 
«  statim, saqout, DETECT(4) ) 

End  if 

Get  channel  calibration  and  delay 

calib  «  chan(5) 
delay  **  chan(6) 

Load  up  previous  filtered  data  to  avoid  start-up  ring 
Filter  the  data  to  be  used  in  the  ST A/ LTA  calculation 

Do  40  i  *  1,75 

f  ilteqC  i)  «*  XX(  i) 

Call  Rfil  (NPTS,BCQEF(0,ICHAN(12)),RAW(76), filter) 

Do  42  i  -  1,75 

XX( i)  -  f ilteq( i+150) 

Set  up  display  times  and  files 
Filter  time  delay  is  0.3  sec  for  Rfil 

ravtim  «*  timei-15 .0+delay 
filtim  ■  timei-15 .0+delay-0 .3 
statim  ■  timei-  7  .5+delay-0 ,3-dtsta 

Define  thresholds 

tsrt  -  TRN0N*LTA 
tend  -  TRN0FF*LTA 
th  *  tsrt 

If  (DETN)  th  -  tend 

Set  up  STA  &  LTA  Storage 
Begin  Loop  over  data 

nsta  ■  1 

Do  1000  filptr  -  l.NPTS.S 
cursta  ■  0.0 

11  -  f ilptr-(S*(NS-l) ) 

12  ■  filptr+(S-l) 
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Do  100  i  ■  il,i2 

cursta  “  cursta+Abs(f ilter( i)) 
nata  ■  nsta+1 
statin  ■  statin  ♦  dtsta 

Shift  detection  flag  array; 

Check  for  STA/LTA  >  threshold; 

Do  110  i  -  1.0P-1 
STA(i)  -  STA( i+1) 

QR(i)  -  QR( i+1) 

STAC  OP)  *  cursta 
big  -  .false. 

OR(QP)  -  0 

If  (cursta. ge.th)  Then 
big  *  .true. 

QR(QP)  -  1 
Endif 
qsun  “  0 
Do  112  i  -  1,QP 

qsun  •  qeum+QR( i) 
qqp  “  .false. 

If  (qsun.ge.Q)  qqp  ■  .true. 

Check  for  detection: 

big  ■  true  if  last  STA  >  threshold 

qqp  ■  true  if  q/q'  STA's  >  threshold 

DETN  ■  true  if  a  detection  is  in  progress 

detsrt  -  .false, 
detend  ■  .false. 

Detection  Start? 

If  (big  .and.  .not. DETN  .and.  qqp)  Then 
detsrt  ■  .true, 
valid  *  .true. 

DETN  ■  .true, 
th  *  tend 

ravptr  •  filptr+75 

Call  Refine  (f ilteq, ravptr ,  STA(OP-l),  LTA,  detptr) 
DETECT(l)  -  filtin  +  (detptr-l)*dt 
DETECT( 2)  -  0.0 
Do  120  i  -  1,0P 

If  (STA(i).gt.DETECT(2))  DETECT(2)  -  STA( i) 
Continue 

DETECTC2)  -  DETECT(2)/Float(S*NS) 

DETECT(3)  -  LTA/ F lost (R*S) 

Check  for  detection  on  data  drop-out. 

iknt  -  0 
drop  “  .false. 

Do  125  i  ■  (detptr-(ISR/4)),detptr+ISR 
If  (RAW(i).ne.O)  Then 
iknt  ■  0 
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Else 

iknt  “  iknt+1 
Endif 

If  ( iknt.ge.ISR/ 2)  drop  -  .true. 

125  Continue 

If  (drop)  Then 

STATUS  -  'Drop-out' 

DETECT ( 4)  -  -999.99 
valid  m  .false. 

Goto  150 

Endif 

c 

c...  Load  up  FFT  buffer 
c 

fftptr  *  detptr-(64/2) 

Do  130  i  *  1,64 

j  ■  fftptr+(i-l) 
xrl( i)  ■  RAW( j) 

130  Continue 
c 

c...  Check  for  Big  spikes 

Call  Spkchk  (64,  xrl ,  IDX,  BRATIO,  BSPK) 

If  (BSPK)  Then 

STATUS  -  'Spike^l ' 
ichan(lO)  *  ichan(10)+l 
valid  ■  .false. 

Endif 

c 

c...  Detrend  &  taper  the  data 
c 

Call  Detrnd  (xrl ,64,0) 

Call  Taper  (xrl ,64,0.1) 
c 

c...  Load  up  complex  array  &  do  FFT 
c 

Do  132  i  *  1 ,64 

132  xcl(i)  *  Cmplx(xrl( i) ,0.0) 

Call  Nlogn  (6,xcl,-1.0) 
c 

c...  Load  up  inverse  transform  buffer 
c...  Remove  low  frequency  components  from  FFT  data 
c...  Do  inverse  transform 
c 

Do  134  i  -  1,  64 

134  xc2(  i)  *  xcl(i) 

xc2(l)  “  (0.0, 0.0) 

Do  140  i  -  2.IFCUT 

xc2(i)  ■  (0.0, 0.0) 
xc2(64— ( i-2) )  “  (0.0, 0.0) 

140  Continue 

Call  Nlogn  (6,xc2,+1.0) 


c 

c  •  •  • 


Pull  out  inverse  xformed  data 


142 


Do  142  i  “  1,64 

xr2(i)  ■  Real(xc2( i)) 

c...  Check  for  Small  spike 
c 

If  (.not.BSPK)  Then 

Call  Spkchk  (64,xr2 .1024.SRATI0,  SSPK) 

If  (SSPK)  Then 

STATUS  -  'Spike;_2' 
ichan(ll)  -  ichan(ll)+l 
valid  ■  .false* 

End  if 

Endif 

c 

c...  Compute  the  amplitude  spectra, 
c 

DETECT(5)  -  0.0 
fmax  “  0.0 
Do  144  i  -  1,33 

xrl(i)  *  Sqrt(Real(xcl( i))**2  +  Aimag(xcl( i))**2) 

If  (xrl( i) .ge.fmax)  Then 
fmax  »  xrl( i) 

DETECT(5)  -  Float( i-l)*df 
Endif 

144  Continue 
c 

c...  Compute  the  spectral  ratio 
c 

sumlo  ■  0.0 
Do  146  i  -  IFL0.IFMID 
146  sumlo  ■  sumlo+xrl( i) 

sumhi  “0.0 

Do  148  i  -  IFMID+1 ,IFHI 
148  sumhi  “  8umhi+xrl(i) 

c 

c...  If  valid  detection,  use  Spectral  Ratio  for  Local/Teleseismic 
c...  discrimination.  ALL  detections  vith  logSR  <  -0.15  are  called 
c...  teleseisms, 
c 

DETECT(4)  -  -999.99 

If  (sumlo. gt. 0.0  .and.  sumhi. gt. 0.0) 

&  DETECT(4)  “  Logl0( sumhi/ sumlo) 

If  (DETECT(4) .It. -0.15)  valid  -  .true. 

If  (valid)  Then 

If  (0ETECT(4) .gt.SPECT)  Then 
STATUS  -  'Local' 
ichan(8)  ■  ichan(8)+l 
Else 

STATUS  “  'Teleseism' 
ichan(9)  ■  ichan(9)+l 
Endif 

Endif 

Endif 

c 

c...  Cont  ~  ie  Detection? 
c 
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If  (DETN  .and.  qqp)  th  K  tend 
c 

c...  Detection  End? 
c 

If  (.not. big  .and.  DETN  .and.  .not. qqp)  Then 
detend  *  .true. 

DETN  -  .false, 
th  K  tsrt 
saqout  «  .false. 

If  ( ST ATUS.eq.' Local'  .or.  STATUS. eq. 'Teleseism')  Then 
saqout  ■  .true. 

Endif 

Call  Report  (CHANlD.statim, STATUS, 

&  DETECT(2) ,DETECT(3) ,DETECT( 1) , saqout, DETECT( 4)) 

Endif 
c 

c...  Update  LTA  every  'R'  STA's; 
c 

130  Continue 

If  (nsta.eq.R)  Then 
If  (DETN)  Then 
If  (valid)  Then 

LTA  »  e2sigd*cursta  +  sigdf*LTA 
Else 

valid  ■  .true. 

Endif 

Else 

LTA  **  eZeig  *curata  +  a igf  *LTA 
Endif 

tsrt  -  TRN0N*LTA 
tend  -  TRN0FF*LTA 
nsta  *  0 
Endif 
c 

c...  Update  start  and  end  detection  thresholds, 
c 

th  *  tsrt 

If  (DETN)  th  -  tend 
1000  Continue 
c 

Return 

End 
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rf  il^.(n,b,x,y) 
int  *n; 

float  b[] ,x[] ,y[ ] ; 

{ 

♦def ine  IRP  7 

•tatic  float  a[IRPl  *  { 

1.0000000e+00, 

0.0000000e+00,  -3.0000000 e+00,  0 .0000000e+00, 

3  .0000000e+00 ,  0 .0000000e+00 ,  -1 .0000000e+00} 


> 


register  int  i,j; 
float  sumax,sumby; 

for  ( i  ■  0  ;  i  <  *n  ;  i++)  { 

•umax  ■  x[ il ; 

for  Cj  -  2  ;  j  <  IRP  ;  j  +-  2) 
•umax  +■  atj]*x[i-j]; 

sumby  ■  bl 1 ] *y I i— 1 ] ; 

for  (j  -  2  ;  j  <  IRP  ;  j++) 

aumby  +-  b[  j ] *y [ i- jl ; 

ylij  *  sumax-aumby; 

> 


♦include  <stdio.h> 
toilet-  () 

{ 

ff lush(stdout) ; 
ff lush(stderr) ; 
return; 
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SAMPLE  OUTPUT 


Unix  DPS  —  Version  80.232 

Detecting  on  7  Sites 
Data  Window:  15  Seconds 
Sampling  Rate:  10  Hertz 
STA  Interval:  0.5  Seconds 
LTA  "  :  1.5 

Detect  on  3  of  3 

LTA  Time  Constant  Threshold 
Normal:  30.0  Secs  2.0 

Detection:  3.0  Secs  1.5 

Filter  Characteristics: 

Order:  3  Type:1.2  -  3.0  Hz  Bandpass 
N  A  Prime(N)  B(N) 

0  1.000000e+  0  1.000000.'+  0 

1  0.  e+  0  -1 .122158e+  0 

2  -3.000000e+  0  1.254692e+  0 

3  0.  e+  0  -8 .377432e-  1 

4  3.000000e+  0  6.884214e-  1 

5  0.  e+  0  -2 .241 16 1 e-  1 

6  -1.000000e+  0  8.427375e-  2 

FFT  Parameters: 

Number  of  Points:  64 

Delta  Time:  0.10000  Sec  Delta  Frequency:  0.15625  Hz 

Spectral  Ratio  Bands: 

Low  7  (  0.94  Hz)  to  13  (  1.88  Hz) 

High  14  (  2.03  Hz)  to  25  (  3.75  Hz) 

log  Ratio  for  Reg ional/Telese ismic  test:  0.20 

Spike  Check  Parameters 

Temporal  Coherency:  1.00  Sec 

"Big"  Ratio:  1.50 

"Small"  Ratio:  2.00 

Low  Cutoff  for  Spikes:  3  (  0.31  Hz) 

Run  will  begin  at  first  record 
Run  will  end  at  last  record 
SAQ  File  is  L07419.saq 
Index  File  is  L07419.gi 
Marker  File  is  L07419.mk 
Gram  File  is  L07419.0.g 

Channel  Information 


Site 

Channel 

Site 

# 

DPS 

Ptr 

Chan 

Ptr 

Chan 

Incr 

Calib 

Delay 

Display 

ALK 

20SHBFAK 

3 

307 

13 

32 

0.10470 

-4.00 

Off 

ALK 

20SHUCAK 

3 

307 

16 

32 

0.08860 

0. 

Off 

ALK 

20SHCNAK 

3 

307 

17 

32 

0.10480 

0. 

Off 

ALK 

20SHNJAK 

3 

307 

18 

32 

0.11570 

0. 

Off 

ALK 

20IHTNAK 

3 

307 

19 

32 

0.21860 

0. 

Off 

PWY 

20SHP999 

4 

710 

13 

28 

0.10000 

-4.00 

Off 

ALK 

20IHATAK 

3 

307 

20 

32 

0.61600 

0. 

Off 

Tape  L07419  starts  at  79/091/06:39:44 
First  data  at  91/06:39:44 
DP  started  at  91/06:39:44 
Record  Number  1 
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Name 

T  ime 

Status 

STA 

LTA 

STA/ LTA 

log  SR 

Max  Freq 

BFAK 

06:39:44.5 

Initial 

0 

7 

0. 

0. 

0. 

UCAK 

06:39:44.5 

Initial 

0 

42 

0. 

0. 

0. 

CNAK 

06:39:44.5 

Initial 

0 

40 

0. 

0. 

0. 

NJAK 

06:39:44.5 

Initial 

0 

28 

0. 

0. 

0. 

TNAK 

06:39:44.5 

Initial 

0 

23 

0. 

0. 

0. 

P999 

06:39:44.5 

Initial 

0 

6 

0. 

0. 

0. 

ATAK 

06:39:44.5 

Initial 

0 

45 

0. 

0. 

0. 

BFAK 

06:39:47.7 

Teleseism 

14 

6 

2.29 

-0.60 

0.78 

1 

P999 

06:39:45.6 

Teleseism 

13 

5 

2.73 

-0.47 

0.63 

2 

P999 

06:40:22.3 

Teleseism 

13 

5 

2.66 

-0.34 

0.63 

3 

N01 

06:41:19 

Data  Gap 

ATAK 

06:41:14.8 

Teleseism 

74 

32 

2.33 

-1.01 

0.16 

4 

N01 

06:41:37 

Data  Gap 

TNAK 

06:41:31.4 

Teleseism 

50 

20 

2.47 

-0.76 

0.31 

5 

N01 

06:44:05 

Data  Gap 

N01 

06:44:25 

Data  Gap 

BFAK 

06:44:24.0 

Teleseism 

16 

7 

2.34 

-0.62 

0.63 

6 

TNAK 

06:44:22.7 

Spike  1 

145 

27 

5.32 

0.06 

0.63 

ATAK 

06:44:32.3 

Teleseism 

87 

30 

2.88 

-0.86 

0.63 

7 

BFAK 

06:45:11.1 

Teleseism 

23 

8 

2.73 

-0.77 

0.94 

8 

NO  I 

06:45:50 

Data  Gap 

N01 

06:46:12 

Data  Gap 

UCAK 

06:45:49.9 

Spike_l 

2866 

149 

19.22 

0.23 

0.31 

P999 

06:45:53.8 

Teleseism 

11 

5 

2.13 

-0.24 

0.78 

9 

NJAK 

06:47:45.9 

Spike  _1 

38 

16 

2.39 

-0.01 

0.47 

TNAK 

06:48:11.3 

Teleseism 

50 

19 

2.60 

-0.68 

0.78 

10 

ATAK 

06:48:11.1 

Teleseism 

72 

28 

2.57 

-0.83 

0.78 

11 

NJAK 

06:49:14.7 

Teleseism 

46 

20 

2.33 

0.03 

0.31 

12 

BFAK 

06:49:17.3 

Teleseism 

25 

9 

2.86 

-0.17 

0.78 

13 

NJAK 

06:49:24.7 

Teleseism 

82 

32 

2.56 

0.07 

0.16 

14 

NJAK 

06:49:30.3 

Teleseism 

99 

37 

2.68 

0.10 

0.16 

15 

TNAK 

06:49:55.3 

Teleseism 

56 

21 

2.68 

-0.86 

0.31 

16 

P999 

06:50:05.9 

Teleseism 

11 

5 

2.12 

-0.16 

0.78 

17 

P999 

06:50:15.3 

Teleseism 

16 

7 

2.52 

-0.37 

0.63 

18 

P999 

06:50:24.7 

Teleseism 

24 

8 

2.96 

-0.68 

0.94 

19 

ATAK 

06:50:34.1 

Teleseism 

63 

30 

2.11 

-0.65 

0.47 

20 

N01 

06:50:52 

Data  Gap 

NJAK 

06:51:16.3 

Teleseism 

59 

21 

2.82 

-0.34 

0.31 

21 

BFAK 

06:51:45.6 

Teleseism 

26 

8 

3.16 

-0.67 

0.63 

22 

CNAK 

06:51:56.4 

Teleseism 

41 

17 

2.47 

-0.28 

0.31 

23 

TNAK 

06:51:56.7 

Teleseism 

55 

21 

2.60 

-0.42 

0.31 

24 

NOl 

06:52:27 

Data  Gap 

N01 

06:53:04 

Data  Gap 

P999 

06:53:03.6 

Teleseism 

11 

5 

2.17 

-0.26 

1.09 

25 
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APPENDIX  V 

PWY  Waveform  and  Detection  Markers 
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